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at  a  distance  of  5  X  from  the  cylinder  and  then  propagated  forward 
to  the  receiver  line.  Diameter,  refractive  index  and  attenuation 
constant  of  the  cylinder  are  6  X,  1.1,  and  4.0  db/cm,  respectively. 

The  attenuation  constant  of  the  medium  is  3.5  db/cm . 56 

2.28  Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  a  cylinder 
illuminated  by  a  point  source  3  X  off  axis  and  10  X  away 

from  the  cylinder,  measured  on  the  receiver  line  at  a  distance  of 
15  X  from  the  cylinder.  Solid  line  represents  the  magnitude  of  the 
field  when  it  is  calculated  directly,  while  the  dashed  line 
corresponds  to  the  case  where  the  field  is  first  calculated  at  a 
distance  of  5  X  from  the  cylinder  and  then  propagated  forward  to 
the  receiver  line.  Diameter,  refractive  index  and  attenuation 
constant  of  the  cylinder  are  6  X,  1.1,  and  4.0  db/cm,  respectively. 

The  attenuation  constant  of  the  medium  is  3.5  db/cm . 57 

2.29  Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  the 
two-cylinder  object  shown  in  Figure  2.15a  measured  on  the  receiver 
line  at  a  distance  of  20  X  from  the  origin.  Solid  line  represents 

the  magnitude  of  the  field  when  it  is  calculated  directly,  while 

the  dashed  line  corresponds  to  the  case  where  the  field  is  first 

calculated  at  a  distance  of  10  X  from  the  origin  and  then 

propagated  forward  to  the  receiver  line.  The  refractive  index  of 

the  cylinders  is  1.1 . 58 
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2.30  Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  the 
two-cylinder  object  shown  in  Figure  ?.15a  measured  on  the  receiver 
line  at  a  distance  of  20  X  from  the  origin.  Solid  line  represents 

the  magnitude  of  the  field  when  it  is  calculated  directly,  while 
the  dashed  line  corresponds  to  the  case  where  the  field  is  first 
calculated  at  a  distance  of  10  X  from  the  origin  and  then 
propagated  forward  to  the  receiver  line.  The  refractive  index  of 
the  cylinders  is  1.5 . 

2.31  Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  the 
two-cylinder  object  shown  in  Figure  2.21a  measured  on  the  receiver 
line  at  a  distance  of  2C  X  from  the  origin.  Solid  line  represents 

the  magnitude  of  the  field  when  it  is  calculated  directly,  while 
the  dashed  line  corresponds  to  the  case  where  the  field  is  first 
calculated  at  a  distance  of  10  X  from  the  origin  and  then 
propagated  forward  to  the  receiver  line.  The  refractive  index  of 
the  cylinders  is  1.1 . 

2.32  Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  the 
two-cylinder  object  shown  in  Figure  2.21a  measured  on  the  receiver 
line  at  a  distance  of  20  X  from  the  origin.  Solid  line  represents 

the  magnitude  of  the  field  when  it  is  calculated  directly,  while 
the  dashed  line  corresponds  to  the  case  where  the  field  is  first 
calculated  at  a  distance  of  10  X  from  the  origin  and  then 
propagated  forward  to  the  receiver  line.  The  refractive  index  of 
the  cylinders  is  1.5 . i 

3.1  Illustration  of  the  Fourier  Diffraction  Projection  Theorem . 

3.2  In  conventional  tomography  the  field  is  measured  on  the  receiver 

line  which  is  perpendicular  to  the  direction  of  the  incident 

plane  wave . ' 

3.3  The  Fourier  transform  of  the  scattered  field  on  the  receiver  line 

gives  the  values  of  the  object's  Fourier  transform  over  a 

semi-circle  in  the  Fourier  domain . ‘ 

3.1  By  rotation  of  the  object  over  360  degrees  the  object's  Fourier 
transform  is  recovered  over  a  set  of  semi-circles  inside  the 
circle  v  2 k  . i 
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3.5  The  four  points  used  for  the  bilinear  interpolation  in  the  frequency 


domain  can  be  assumed  to  be  the  edges  of  a  parallelogram  if  the 
rotational  angular  step  is  small  and  the  number  of  points  on  the 
receiver  line  is  large . 79 

3.6  A  schematic  of  the  geometry  of  the  Synthetic  Aperture  Diffraction 

Tomography  system . 80 

3.7  The  direction  of  the  wave  vectors  K  and  A  relative  to  the 

space  coordinates  are  depicted . 82 

3.8  For  a  fixed  orientation  of  the  object  in  the  transmission  mode,  the 

object's  Fourier  transform  is  obtained  in  the  region  shown . 85 


3.9  For  two  rotational  positions  of  the  object,  90  degrees  apart,  the 
object's  Fourier  transform  is  recovered  inside  a  circle  of  radius 


v2k0.  The  coverage  in  the  region  C  is  double . 86 

3.10  The  size  of  the  object,  S,  the  distance  of  the  object  from  the 
transmitter  line,  d,  and  the  distance  between  the  arrays,  D, 
are  depicted  in  this  figure . 88 


3.11  When  the  transmitter  array  in  figure  3.6  is  utilized  only  in  the 

reflection  mode,  for  a  fixed  orientation  of  the  object,  the  object's 
Fourier  transform  is  recovered  inside  the  region  shown . 90 

3.12  The  coverage  inside  the  circle  of  radius  2k0  is  achieved  by  utilizing 

the  system  both  in  the  transmission  and  the  reflection  modes  for 
only  one  rotational  position  of  the  object . 91 

3.13  In  the  reflection  mode,  the  object’s  Fourier  transform  is  obtained 

over  the  semi-circles  SC(— kx ky ),  one  of  which  is  the  semi-circle 
centered  at  point  A.  By  moving  this  point  on  the  dashed 
semi-circle,  the  coverage  shown  in  figure  3.11  is  achieved . 93 

3.1  1  The  points  A  and  B  are  the  center  of  semi-circles  SC'(-kx,-ky) 
which  pass  over  the  point  D.  Point  D  is  restricted  to  the  region 
shown . 96 
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3.15  The  curved  surface  shows  the  coven  go  achieved  in  the  object's 

“Fourier-Laplace"  domain  when  the  medium  is  attenuatin  g.  The 
attenuation  of  the  medium  was  taken  to  be  3.0  db/cm . 100 

3.16  For  a  lossy  medium,  a  selective  part  of  the  collected  data  n  the 

Synthetic  Aperture  technique,  and  many  rotational  views  of  the 

object,  result  in  the  values  of  the  object's  Fourier  transform 

over  radial  lines . 102 

4.1  Cross  section  of  a  two-component  object  reconstructed  using  the 

conventional  diffraction  reconstruction  technique.  Diameter  of 

the  cylinders  is  6  X  and  their  refractive  index  is  1.02. 

a)  Only  the  first-order  scattered  fields  are  used  for  generating 

the  data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 

included  in  the  projection  data  for  this  reconstruction . 108 

4.2  Cross  section  of  a  two-component  object  reconstructed  using  the 

conventional  diffraction  reconstruction  technique.  Diameter  of 

the  cylinders  is  6  X  and  their  refractive  index  is  1.03. 

a)  Only  the  first-order  scattered  fields  are  used  for  generating 

the  data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 

included  in  the  projection  data  for  this  reconstruction . 109 

4.3  Cross  section  of  a  two-component  object  reconstructed  using  the 

conventional  diffraction  reconstruction  technique.  Diameter  of 

the  cylinders  is  6  X  and  their  refractive  index  is  1.05. 

a)  Only  the  first-order  scattered  fields  are  used  for  generating 

the  data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 

included  in  the  projection  data  for  this  reconstruction . 110 

4.4  Cross  section  of  a  three-component  object  reconstructed  using  the 

conventional  diffraction  reconstruction  technique.  Diameter  of 

the  cylinders  is  6  X  and  their  refractive  index  is  1.05. 

a)  Only  the  first-order  scattered  fields  are  used  for  generating 

the  data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 

included  in  the  projection  data  for  this  reconstruction . Ill 
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(a)  The  condition  when  the  object  components  block  each  other  is 
illustrated  here,  (b)  With  the  oblique  reception  illustrated  here 
the  “direct  blockage”  is  avoided.  Direct  blockage  occurs  when 
the  direction  of  the  propagation  of  the  incident  field  is  lined 
up  with  the  line  joining  the  two  object  components  and  the 
normal  to  the  receiver  array . 113 


4.6  In  the  second  approach  several  plane  waves  are  used  to  illuminate 
the  object.  The  frequency  domain  coverage  obtained  with  the  three 
plane  waves  shown  here  are  averaged  prior  to  inversion . 114 


4.7  (a)  When  the  incident  plane  wave  is  normal  to  the  receiver  line, 

the  Fourier  transform  of  the  projection  gives  the  values  of  the 
object’s  Fourier  transform  over  a  semi-circle  as  shown  here. 

(b)  When  the  incident  wave  is  at  an  oblique  angle  with  the 
receiver  line,  the  coverage  is  still  on  a  semi-circle,  although 
positioned  differently . 116 

4.8  In  the  third  approach,  again  a  number  of  plane  waves  are  used  to 

illuminate  the  object,  however  no  independent  frequency  domain 
coverages  are  generated  from  each  of  these  plane  waves.  For  each 
rotational  position  of  the  object,  the  coverage  generated  by  all 
the  illuminating  waves  is  extracted  first,  and  then  the  object 
is  rotated  till  the  frequency  domain  is  filled  up . 119 

4.9  Reconstructions  of  a  two-component  object  using  the  first  approach. 

The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and 
positioned  at  a  distance  of  12  X  from  each  other.  The  projections 
are  sampled  at  64  points  with  a  sampling  interval  of  0.45  X.  The 
receiver  line  is  20  X  from  the  center  of  the  coordinates.  The  angle 
between  the  direction  of  illumination  and  the  normal  to  the  receiver 
line  is:  (a)  0°;  (b)  7.5°;  (c)  15°;  (d)  30° . 120 
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4.10  Reconstructions  of  a  two-component  object  using  the  second 

approach.  The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter, 
and  positioned  at  a  distance  of  12  X  from  each  other.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of 
0.45  X.  The  receiver  line  is  20  X  from  the  center  of  the  coordinates. 

(a)  3  incident  plane  waves,  with  the  angle  of  illumination  between 

-15°  and  15°.  (b)  3  incident  plane  waves,  with  the  angle  of 
illumination  between  -30°  and  30° . 123 

4.11  Reconstructions  of  a  two-component  object  using  the  third 

approach.  The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter, 
and  positioned  at  a  distance  of  12  X  from  each  other.  The  projections 
are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X.  The 
receiver  line  is  20  X  from  the  center  of  the  coordinates,  (a)  3  incident 
plane  waves,  with  the  angle  of  illumination  between  -7.5°  and  7.5° 
and  32  receiver  positions,  (b)  5  incident  plane  waves,  with  the 
angle  of  illumination  between  -15°  and  15°  and  32  receiver  positions. 

(c)  5  incident  plane  waves,  with  the  angle  of  illumination  between 
-30°  and  30°  and  16  receiver  positions . 124 

4.12  Reconstructions  of  a  two-component  object  using  the  first  approach 

consisting  of  oblique  reception.  The  cylinders  of  refractive  index 
1.08  are  4  X  in  diameter,  and  positioned  at  a  distance  of  8  X  from 
each  other.  The  projections  are  sampled  at  64  points  over  a 
sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from  the 
center  of  the  coordinates.  The  angle  between  the  direction  of 
illumination  and  the  normal  to  the  receiver  line  is:  (a)  0°; 

(b)  10°;  (c)  20°;  (d)  30° . 126 

4.13  Reconstructions  of  a  two-component  object  using  the  second 

approach.  The  cylinders  of  refractive  index  1.08  are  4  X  in  diameter, 
and  positioned  at  a  distance  of  8  X  from  each  other.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of 
0.45  X.  The  receiver  line  is  20  X  from  the  center  of  the  coordinates. 

(a)  3  incident  plane  waves,  with  the  angle  of  illumination  between 
-10°  and  10°.  (b)  5  incident  plane  waves,  with  the  angle  of 
illumination  between  -20°  and  20°.  (c)  7  incident  plane  waves, 
with  the  angle  of  illumination  between  -30°  and  30° 


128 


4.14  Reconstructions  of  a  two-component  object  using  the  third 

approach.  The  cylinders  of  refractive  index  1.08  are  4  X  in 
diameter,  and  positioned  at  a  distance  of  8  X  from  each  other. 

The  projections  are  sampled  at  64  points  over  a  sampling  interval 
of  0.45  X.  The  receiver  line  is  20  X  from  the  center  of  the 
coordinates,  (a)  5  incident  plane  waves,  with  the  angle  of 
illumination  between  -20°  and  20°  and  16  receiver  positions. 

(b)  7  incident  plane  waves,  with  the  angle  of  illumination 

between  -30°  and  30°  and  8  receiver  positions . 130 

4.15  Reconstructions  of  a  three-component  object  using  the  first 

approach  consisting  of  oblique  reception.  The  cylinders  of 
refractive  index  1.05  are  4  X  in  diameter,  and  positioned  at 
a  distance  of  6  X  from  the  center  of  coordinates.  The  projections 
are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X.  The 
receiver  line  is  20  X  from  the  center  of  the  coordinates.  The 
angle  between  the  direction  of  illumination  and  the  normal  to  the 
receiver  line  is  (a)  0°,  (b)  10°,  (c)  20°,  (d)  30° . 132 

4.16  Reconstructions  of  a  three-component  object  using  the  second 

approach.  The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter, 
and  positioned  at  a  distance  of  6  X  from  the  center  of  coordinates. 

The  projections  are  sampled  at  64  points  over  a  sampling  interval 
of  0.45  X.  The  receiver  line  is  20  X  from  the  center  of  the 
coordinates,  (a)  3  incident  plane  waves,  with  the  angle  of 
illumination  between  -10°  and  10°.  (b)  5  incident  plane  waves, 
with  the  angle  of  illumination  between  -20°  and  20°.  (c)  7  incident 
plane  waves,  with  the  angle  of  illumination  between  -30°  and  30°.  ...134 

4.17  Reconstructions  of  a  three-component  object  using  the  third 

approach.  The  cylinders  of  refractive  index  1.05  are  4  X  in 
diameter,  and  positioned  at  a  distance  of  6  X  from  the  center 
of  coordinates.  The  projections  are  sampled  at  64  points  over  a 
sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from  the 
center  of  the  coordinates,  (a)  5  incident  plane  waves,  with  the 
angle  of  illumination  between  -20°  and  20°  and  16  receiver 
positions,  (b)  7  incident  plane  waves,  with  the  angle  of  illumination 
between  -3 0°  and  30°  and  16  receiver  positions . 136 
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4.18  A  simple  interpolation  scheme  for  the  Synthetic  Aperture  technique. 

Points  A  and  A' ,  B  and  B' ,  and  C  and  C  lying  on  two  different 
semi-circles  SC(-kx,-ky)  posses  equal  values  of  respectively. 

(a)  small  values  of  ky.  (b)  values  of  ky  close  to  k0 . 138 

4.19  Points  Ci  and  C2  are  the  intersection  points  of  the  semi-circles 

04  and  02  with  the  circle  centered  at  the  origin  going  through 
point  C.  The  neighboring  points  of  Cj  and  C2  are  chosen  for  the 
kind  of  interpolation  desired . 140 

4.20  Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with 

diameter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic 
Aperture  technique.  Transmitter  and  receiver  arrays,  positioned 
20  X  apart,  have  32  elements  with  inter-element  distance  of  0.45  X. 

(a)  Simple  nearest  neighbor  interpolation,  (b)  Modified  nearest 
neighbor  interpolation,  (c)  Simple  bilinear  interpolation. 

(d)  Modified  bilinear  interpolation . 141 

4.21  Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with 

diameter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic 
Aperture  technique.  Transmitter  and  receiver  arrays,  positioned 
20  X  apart,  have  64  elements  with  inter-element  distance  of  0.45  X. 

(a)  Simple  nearest  neighbor  interpolation,  (b)  Modified  nearest 
neighbor  interpolation,  (c)  Simple  bilinear  interpolation. 

(d)  Modified  bilinear  interpolation . 143 

4.22  Magnitude  of  the  scattered  field  from  a  lossless  cylinder,  embedded 

in  a  lossless  medium,  as  a  function  of  transmitter  and  receiver 

positions.  Diameter  and  refractive  index  of  the  cylinder  are  6  X 

and  1.03,  respectively . 146 

4.23  Magnitude  of  the  Fourier  transform  of  the  scattered  field  from  a 

lossless  cylinder,  embedded  in  a  lossless  medium,  with  respect 
to  transmitter  and  receiver  position.  Diameter  and  refractive 
index  of  the  cylinder  are  6  X  and  1.03,  respectively,  (a)  32  element 
transmitter  and  receiver  arrays,  (b)  128  element  transmitter  and 
receiver  arrays . 147 
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4.24  Due  to  the  loss  of  high  frequency  components  in  the  projections  in 

the  Synthetic  Aperture  technique,  for  one  rotational  view  of  the 

object  the  object’s  Fourier  transform  is  obtained  only  in  the 

region  shown . 148 

4.25  When  the  high  frequency  information  in  the  projections  is  lost, 

larger  number  of  views  can  increase  the  resolution  and  decrease 
the  distortion  in  the  reconstructed  image.  The  coverage  shown 
corresponds  to  three  rotational  views  of  the  object  60°  apart . 149 

4.26  Magnitude  of  the  reconstructed  object’s  Fourier  transform  in  a 

lossless  medium.  The  object  is  a  lossless  cylinder  with  diameter 
and  refractive  index  of  6  X  and  1.03,  respectively,  (a)  32  element 
transmitter  and  receiver  arrays,  (b)  128  element  transmitter  and 
receiver  arrays . 151 

4.27  Real  part  of  the  reconstructed  object  function  of  a  cylinder  with 

diameter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic  - 

Aperture  technique.  Transmitter  and  receiver  arrays,  with  the 
inter-element  distance  of  0.45  X,  are  positioned  20  X  apart,  (a)  32 
element  arrays,  (b)  64  element  arrays,  (c)  128  element  arrays . 152 

4.28  Imaginary  part  of  the  reconstructed  object  function  (difference 

in  attenuation  constant  of  object  and  medium)  of  a  cylinder  with 
diameter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic 
Aperture  technique.  Transmitter  and  receiver  arrays,  with  the 
inter-element  distance  of  0.45  X,  are  positioned  20  X  apart,  (a)  32 
element  arrays,  (b)  64  element  arrays,  (c)  128  element  arrays . 154 

4.29  Magnitude  of  the  reconstructed  refractive  index  of  a 

two-component  object,  using  the  Synthetic  Aperture  technique. 

The  components  are  cylinders  with  diameter  of  6  X  and  refractive 
index  of  1.03.  Transmitter  and  receiver  arrays  positioned  30  X 
apart,  have  128  elements  with  the  inter-element  distance  of  0.45  X. 
a)  Only  the  first-order  scattered  fields  are  used  for  generating  the 
data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 
included  in  the  projection  data  for  this  reconstruction . 156 
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4.30  Magnitude  of  the  reconstructed  refractive  index  of  a 

two-component  object,  using  the  Synthetic  Aperture  technique. 

The  components  are  cylinders  with  diameter  of  6  X  and  refractive 
index  of  1.05.  Transmitter  and  receiver  arrays  positioned  30  X 
apart,  have  128  elements  with  the  inter-element  distance  of  0.45  X. 
a)  Only  the  first-order  scattered  fields  are  used  for  generating  the 
data  for  this  reconstruction,  b)  Doubly  scattered  fields  are 
included  in  the  projection  data  for  this  reconstruction . 157 

4.31  Magnitude  of  the  reconstructed  refractive  index  of  a  three-component 

object,  using  the  Synthetic  Aperture  technique.  The  components 

are  cylinders  with  diameter  of  6  X  and  refractive  index  of  1.05. 

Transmitter  and  receiver  arrays  positioned  30  X  apart,  have 

128  elements  with  the  inter-element  distance  of  0.45  X.  a)  Only 

the  first-order  scattered  fields  are  used  for  generating  the  data 

for  this  reconstruction,  b)  Doubly  scattered  fields  are  included  in 

the  projection  data  for  this  reconstruction . 158 

4.32  Energy  loss  of  plane  wave  components  resulting  from  an  angular 
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4.35  Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive 

index  of  1.01,  using  the  conventional  technique.  Both  the  object  and 
the  medium  are  lossless.  The  projections  are  sampled  at  64  points 
over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  at  a  distance 
of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference  in 
attenuation  constant  of  the  object  and  the  medium . 166 
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0.45  X.  The  receiver  line  is  at  a  distance  ol  10  X  from  the  cylinder. 
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4.40  Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive 

index  of  1.03,  using  the  conventional  technique.  The  attenuation 
constant  is  3.0  db/cm  in  both  the  object  and  the  medium.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of 
0.45  X.  The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder. 

(a)  Refractive  index,  (b)  Difference  in  attenuation  constant  of 

the  object  and  the  medium . 171 
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of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference 
in  attenuation  constant  of  the  object  and  the  medium . 172 

4.42  Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive 

index  of  1.0,  using  the  conventional  technique.  The  attenuation  : 

constant  of  the  medium  and  the  object  are  2.0  db/cm  and 
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4.45  Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive 

index  of  1.01,  using  the  conventional  technique.  The  medium  is 
lossless  while  the  attenuation  constant  of  the  object  is  0.5  db/cm. 

The  projections  are  sampled  at  64  points  over  a  sampling 

interval  of  0.45  X.  The  receiver  line  is  at  a  distance  of  10  X 

from  the  cylinder,  (a)  Refractive  index,  (b)  Difference  in 

attenuation  constant  of  the  object  and  the  medium . 177 
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4.53  Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive 
index  of  1.02,  using  the  Synthetic  Aperture  technique.  The 
attenuation  constants  of  the  object  and  medium  are  3.5  db/cm 
and  3.0  db/cm,  respectively.  The  transmitter  and  the  receiver 
arrays,  positioned  20  X  apart,  have  128  elements  with  the 
inter-element  distance  of  0.15  X.  The  attenuation  of  the 
medium  has  been  compensated  for.  (a)  refractive  index  of  the 
object  (b)  difference  in  attenuation  constant  of  the  object 
and  the  medium . 
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CHAPTER  1 
INTRODUCTION 


Tomography  refers  to  the  measurement  in  cross-sect o  planes  uf  t:-> 
parameters  from  either  transmission  or  reflection  data.  I  i  most  ce_ -s,  i 
?is:-n«  structure  i-  illuminated  from  many  different  direct  i.  n-  and  the  image 
two-dimensional  map  of  the  chosen  ti.-.-ue  parameter)  is  rec  a  .■meted  from  I 
data  •,  ollected  either  in  transmission  or  reflection  (called  project:  ,ns). 

Tfmre  is  a  fundamental  diiTerence  between  tomographic  imaging  with 
rays  on  the  one  hand,  and  with  ultrasound  and  microwaves  <n  the  other  hat 
X-ravs.  being  non-diffracting,  travel  in  straight  lines,  and  therefore,  t 
projection  data  measure  the  line  integral  of  some  object  parameter  aln 
straight  lines.  This  makes  it  possible  to  apply  the  Fourier  Slice  '!  h»or< 
[Rosen felds-jj,  which  say  %  that  the  Fourier  transform  of  a  pr>T  ti  :  :s  e.jual 
a  slice  of  the  two-dimensional  Fourier  transform  of  the  objt  •••.  tin- 
f'  rms  the  foundation  of  the  filtered-backproj-ctjon  algoro  m.;  ■  ;  with  r: 

diffracting  sources,  such  as  x-rays. 

On  the  other  hand,  when  either  microwaves  or  ultrasound  are  used  I 
tomographic  imaging,  the  energy  often  does  not  propagate  along  straight  1m 
When  the  object  inhomogeneities  arc  large  compared  to  a  v,  iv  .  length,  cmt 
propagation  is  characterized  by  refraction  and  multipath  •  ■<  N * ■ " r . 

amounts  of  ray  bending  induced  by  refraction  can  be  taken  uo  account 
combining  algebraic  reconstruction  algorithms  [AnderscnMj  w;;h  iigital  r 
tracing  and  ray  linking  algorithms  (Andersen^'ij.  The  multipath  effects,  can-' 
by  the  arrival  at  the  receiving  transducer  of  more  than  ou<  ray  prnpagati 
through  different  parts  of  the  refracting  object,  can  usually  !•’  elmiinat 
through  homomorphic  processing  of  the  projection  data  [(  rawArds. 

W’mn  t  J,i-  object  inhomogeneities  become  comparable  m  -mo  t<> 
wave  eiigth,  it  is  not  even  appropriate  to  talk  about  propagation  along  l.n 
and  rays,  and  energy  transmission  must  be  discussed  in  terms  of  wa.-  fr  a 
and  fo  lds  scfitt.-rcd  bv  the  inhomogeneities.  With  diffracting  illurnmat  i  u: 
e\e.  t  -oitition  to  the  problem  is  made  difficult  by  the  complicated  nature  of  * 
’vv  itter  interaction.  These  interactions  are  descri*  d  by  ditfereip ;  i! 
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integral  equations  which,  in  general,  <lu  not  possess  closed-form  soim Two 
types  of  solutions,  called  the  first-order  Born  and  Rytov,  are  avail. a!  i-  :.t  this 
time  if  one  can  get  away  with  the  assumption  of  weak  interactions  implying 
small  inhomogeneities.  By  either  theoretical  reasoning  or  computer  simulations 
on  simple  cross-sectional  distributions,  it  is  usually  possible  to  determine  the 
magnitude  of  the  inhomogeneities,  which  if  exceeded  would  lead  to  a 
breakdown  in  the  first-order  Born  or  Rytov  based  solutions. 

Small  perturbation  (Born  and  Rytov)  solutions  for  reconstruct'!!  n  were 
initiated  in  1960  by  Wolf  [WolfSOj.  He  showed  that  with  an  incident 
monochromatic  plane  wave,  the  Fourier  transform  of  the  scattered  data 
measured  on  a  receiving  planar  aperture  and  the  Fourier  transform  of  the 
object,  function  are  linearly  related.  Similar  results  for  three-dimensional 
reconstruction  were  obtained  by  Dandliker  and  Weiss  [Dandliker'/O],  Carter 
used  this  method  for  reconstructing  refractive  index  distribution  by  using 
holograms  to  record  the  scattered  fields  in  the  Fresnel  zone  of  a 
semitransparent  object  [CarterTO],  In  spite  of  many  experimental  problems,  his 
results  were  quite  promising.  Iwata  and  Nagata  addressed  some  of  the 
limitations  of  this  approach  and  proposed  a  cylindrical  receiver  line  vis-a-vis 
the  straight  line  (as  suggested  by  Wolf)  for  small  objects  [Iwata75], 

Some  theoretical  procedures  for  reconstructing  weakly  inhomogeneous 
objects  hate  been  recently'  derived  by  several  researchers.  A  three-dimensional 
inverse  scattering  scheme  for  a  spherical  recording  aperture  enclosing  the  object 
was  derived  by  Ball  and  Stenger  [BallSO] .  They  considered  monochromatic 
spherical  wave  sources  for  the  object  insonification  and  point  receivers  for  the 
measurement  of  the  scattered  field.  Norton  and  Linzer  [NortouSl]  have 
considered  broadband  insonification  to  derive  reconstruction  methods  for 
planar,  spherical  and  cylindrical  geometries. 

The  method  suggested  by  Wolf  has  been  employed  by  Mueller  et  ai. 
[Mueller?*).  KavehTOaj  to  implement  a  cross-sectional  reconstruction  scheme  for 
imaging  the  refractive  index  of  objects.  They  have  shown  that  a  coverage  of 
the  Fourier  transform  of  the  object  function  can  be  achieved  by  measuring  the 
scattered  field  on  a  receiving  line  for  different  rotational  positions  of  the  object 
'iver  -%0  degrees.  Mueller  et  al.’s  method  gives  the  values  of  the  Fourier 
tran-form  of  the  object  on  a  set  of  arcs,  while  for  reconstructing  the  object 
funeti'  n  through  the  discrete  Fourier  transform  operation  on  a  digital 
compiler.  the  values  are  required  on  a  two-dimensional  rectangular  grid. 
J)evnney  if),  vaneyvjj  has  derived  a  method  (suggested  independently  by 
Soumekh  e-,  -j.  [Soumrkh^.'tJ  1,  which  he  has  called  the  Fourier  back- 


propagation  technique  to  solve  this  problem.  This  approach  i-  similar  in  spirit 
to  Fourier  Back-Projection  for  straight  ray  tomography.  A  drawback  of  this 
approach  is  its  extensive  demand  on  computer  time. 

Nahamoo  and  Kak  [Nahamoo82]  generalized  the  technique  of  diffraction 
tomography  as  presented  in  [Mueller78j  to  derive  a  new  scheme  which,  in 
principle,  permits  the  use  of  any  type  of  insonification,  as  opp  i  to  only  the 
plane  wave  type  which  is  experimentally  difficult  to  generate.  With  this 
approach,  the  requirement  that  the  object  be  viewed  over  360  degrees  is 
replaced  by  the  need  for  merely  two  views,  which  ideally  should  be  90  degrees 
apart  (however,  that  is  not  strictly  necessary).  For  each  of  these  two  views  the 
method  requires  independent  traversal  movements  for  both  the  transmitter  and 
the  receiver  on  two  parallel  lines  in  the  plane  of  the  cross-section. 

Kenue  and  Greenleaf  [Kenue82]  extended  the  approach  of  Mueller  et  al. 
[Mucller79]  to  the  multi-frequency  case.  They  have  shown  that  by  measuring 
the  projections  at  more  than  one  frequency  one  can  reduce  the  number  of  views 
required  and  thereby  higher  speeds  in  data  collection  may  be  achieved.  Note 
that  the  total  amount  of  data  collected  remains  basically  the  same,  because  for 
each  rotational  position  of  the  object,  each  projection  must  now  be  recorded  at 
different  frequencies.  It  is  possible  that  in  the  actual  implementation,  the 
frequency  dependence  of  tissue  attenuation  may  become  a  source  of  difficulty 
for  this  method. 

Recently  Johnson  et  al.  [Johnson83,  Tracy83,  Johnson84]  have  -uggested  a 
method  based  on  solving  systems  of  nonlinear  algebraic  equations  that  are 
derived  by  applying  the  method  of  moments  to  a  Sine  basis  function  expansion 
of  the  fields  and  scattering  potential.  Although  it  seems  that  the  limitation  of 
this  approach  is  quite  similar  to  that  of  the  Born  approximation,  the  domain  of 
applicability  of  this  approach  still  remains  to  be  proved. 

Diffraction  tomography  has  the  potential  of  providing  us  with  techniques 
for  the  quantitative  measurement  of  tissue  parameters  that  take  into  account 
the  inherently  diffractive  nature  of  wave  propagation.  Diffraction  tomography 
methodology  is  not  limited  by  the  assumptions  of  geometric  propagation,  since 
concepts  such  as  rays  are  not  invoked,  and  is  capable  of  providing  us  with 
distributions  for  virtually  all  the  tissue  parameters.  Unfortunately,  at  this  time 
the  promise  of  this  approach  can  only  be  realized  for  tissues  that  are  weakly 
scattering;  but  it  is  important  to  bear  in  mind  that  the  condition  of  weak 
.-cat t ering  is  not  a  theoretical  limitation,  and  work  is  progressing  at  many 
ro'e'.urch  centers  toward  the  development  of  computational  algorithms  for 
handling  larger  tissue  inhomogeneities. 
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Chapter  2  is  devoted  to  the  forward  scattering  process.  Propagation  of 
the  angular  spectrum  of  a  field  in  a  homogeneous  medium  is  presented  in 
section  2.1.  In  section  2.2  scattering  from  single  component  objects  is 
discussed.  The  scattering  process  is  considered  for  the  case  of  plane  wave 
incidence  as  well  as  point  source  illumination.  The  role  of  attenuation  in  the 
object  as  well  as  in  the  medium  is  described  and  simulation  results,  based  on 
the  theoretical  formulation  presented,  are  shown.  In  section  2  3,  a  new 
procedure  for  calculation  of  scattered  fields  from  a  multi-component  object  is 
proposed  and  computer  simulation  results  presented  and  discussed.  The 
prospect  of  the  new  procedure  in  an  attenuation  medium  is  examined  at  the 
end  of  that  section.  Finally  in  section  2.4,  a  testing  procedure  for  the 
calculated  scattered  fields  is  used  to  check  the  results  of  the  previous  sections. 

In  chapter  3  the  inverse  scattering  problem  is  considered.  The  formulation 
of  the  process  is  presented  in  section  3.1,  and  the  two  well  known  Bom  and 
Rytov  approximations  discussed.  The  conventional  Diffraction  Tomography 
and  related  reconstruction  techniques  are  discussed  in  section  3.2.  The 
derivation  of  the  Synthetic  Aperture  Technique  in  the  transmission  mode  is 
reviewed  in  section  3.3,  and  the  corresponding  theory  for  the  reflection  mode  as 
well  as  the  coverage  achieved  in  a  lossy  medium  are  derived. 

A  simulation  study  of  multiple  scattering  and  attenuation  phenomena  in 
the  inverse  problem  is  presented  in  chapter  4.  In  section  4.1,  artifacts  caused 
by  multiple  scattering  are  studied  by  considering  computer  simulation  results. 
Several  approaches  to  reduce  these  artifacts  are  also  discussed.  A  simulation 
study  of  attenuation  in  Diffraction  Tomography  is  presented  in  section  4.2. 

Although  the  only  source  of  energy  we  have  considered  in  our  study  is 
ultrasound,  the  formulation  of  the  problem  is  the  same  for  microwave  energy. 
The  pressure  field  has  to  be  replaced  by  the  scalar  potential  in  the  case  of 
microwave.  However,  the  incident  electric  field  has  to  be  polarized  either  along 
or  perpendicular  to  the  axis  of  the  cylindrical  objects  considered.  Any  further 
extension  of  the  theory  still  remains  to  be  explored. 
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CHAPTER  2 
FORWARD  PROBLEM 


When  a  wave  encounters  an  obstacle,  some  part  of  the  wave  energy  is 
deflected  from  its  original  course.  It  is  usual  to  define  the  difference  between 
the  actual  wave  and  the  undisturbed  wave,  which  would  be  present  if  the 
obstacle  was  not  there,  as  the  scattered  wave.  When  a  plane  wave,  for 
instance,  strikes  a  body  in  its  path,  in  addition  to  the  undisturbed  plane  wave 
there  is  a  scattered  wave,  spreading  out  from  the  obstacle  in  all  directions, 
distorting  and  interfering  with  the  plane  wave.  If  the  obstacle  is  very  large 
compared  with  the  wavelength,  half  of  this  scattered  wave  spreads  out  more  or 
less  uniformly  in  all  directions  from  the  scatterer,  and  the  other  half  is 
concentrated  behind  the  obstacle  in  such  a  manner  as  to  interfere  destructively 
with  the  unchanged  plane  wave  behind  the  obstacle,  creating  a  sharp-edged 
shadow  there.  This  is  the  case  of  geometric  optics;  in  this  case  the  half  of  the 
scattered  wave  spreading  out  uniformly  is  called  the  reflected  wave,  and  the 
half  responsible  for  the  shadow  is  called  the  interfering  wave  [Morse68],  If  the 
obstacle  is  very  small  compared  with  the  wavelength,  the  scattered  wave  is 
propagated  out  in  all  directions,  and  there  exists  no  sharp-edged  shadow.  In 
the  intermediate  cases,  where  the  obstacle  is  about  the  same  size  as  the 
wavelength,  a  variety  of  curious  interference  phenomena  can  occur. 

In  our  studies,  we  have  treated  the  forward  scattering  process,  as  weii  as 
the  inverse  process,  as  a  two-dimensional  problem.  In  other  words,  the  objects 
under  consideration  are  assumed  to  be  invariant  in  one  dimension.  Although 
most  of  the  Diffraction  Tomography  formulations  and  reconstruction  techniques 
are  or  can  be  generalized  to  the  three-dimensional  case,  we  believe  that  most  of 
the  current  effort  should  be  directed  toward  resolving  the  many  important 
limitations  still  facing  the  two-dimensional  problem. 

If  the  forward  problem  can  be  represented  as  a  boundary  value  problem,  it 
is  possible  to  find  an  exact  solution  which  consists  of  writing  down  a  coordinate 
system  matching  the  boundaries  of  the  system,  and  then  using  the 
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eigenfunction  expansions  for  solving  the  partial  differential  equations. 
Unfortunately,  only  the  circular  inhomogeneities  can  be  solved  for  in  this 
manner.  In  the  past,  this  procedure  has  been  the  only  tool  for  the  generation 
of  simulated  data  for  diffraction  tomography  [Iwata75,  Kaveh79b,  Azimi83], 
But  note  the  limitation:  the  method  is  suitable  only  for  single  component 
objects  that  are  circular  in  shape. 

However,  if  the  object  inhomogepeities  are  very  small  such  that  the  well 
known  Born  or  Rytov  approximations  could  be  justified,  the  scattered  held  can 
be  calculated  by  modeling  the  scattering  process  as  a  linear  system  [Slanev84]. 
Use  of  iterative  techniques,  such  as  higher  order  Born  or  Rytov  techniques,  can 
increase  the  domain  of  applicability  of  this  approach  to  some  extent. 

Since  the  Angular  Spectrum  expansion  has  been  extensively  utilized  in  our 
study,  either  in  the  calculation  of  scattered  fields  or  in  the  development  of 
reconstruction  techniques,  section  2.1  is  devoted  to  the  propagation  of  the 
angular  spectrum  of  a  field  in  either  lossy  or  lossless  medium.  Section  2.2  is 
devoted  to  the  solution  of  the  scattered  field  from  a  single  object  illuminated 
either  by  a  plane  wave  or  a  point  source.  Related  computer  simulation  results 
are  also  presented  in  that  section. 

When  an  object  possesses  more  than  one  component,  due  to  the 
interaction  between  the  components  the  totai  scattered  fields  cannot  be 
considered  as  a  superposition  of  the  fields  due  to  each  of  the  components 
individually,  as  generated  by  the  above  method.  In  section  2.3,  we  present  a 
new  way  of  including  the  interactions  between  the  components  by  expanding 
the  scattered  fields  from  each  component  as  a  superposition  of  plane  waves  and 
then  analyzing  the  interaction  of  each  one  of  these  plane  waves  with  the  other 
object  components,  and  finally  synthesizing  the  total  field.  The  effect  of 
component  interaction  on  the  scattered  fields  is  discussed  in  this  section  by- 
considering  the  numerous  computer  simulation  results  presented.  Finally,  at 
the  end  of  this  section,  the  prospect  of  this  new  technique  in  an  attenuating 
medium  is  examined.  In  section  2.1  we  have  considered  a  testing  procedure  for 
checking  the  results  of  our  computer  programs. 

Past  efforts  on  the  forward  problem  as  applied  to  Diffraction  Tomography 
have  ignored  the  attenuation  present  in  the  medium  or  the  object.  Although 
ultrasound  undergoes  only  a  slight  attenuation.,  in  microwave  imaging 
attenuation  in  the  medium  and  the  object  can  be  a  limiting  factor.  :  this 
chapter,  v- 1 •  have  also  talked  about  the  extension  of  the  algorithms  for  the 
calculation  of  scattrrv'  fields  when  both  the  object  and  the  m  bum  are  f.ssv. 
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2.1  Propagation  of  Angular  Spectrum 

Formulating  the  propagation  of  fields  in  the  Fourier  domain  is  a  useful 
concept.  This  is  especially  true  from  a  computational  point  of  view,  since  FFT 
is  computationally  very  efficient.  In  a  lossless  medium  the  field  and  its  plane 
wave  spectrum  are  related  through  a  Fourier  transform,  thus  each  plane  wave 
can  be  propagated  to  a  new  receiver  line  by  the  addition  of  a  phase  shift.  On 
the  other  hand,  in  a  lossy  material  the  Fourier  transform  does  not  represent  the 
field  in  a  form  that  is  physically  realistic  and  we  must  follow  a  pure 
mathematical  representation  of  the  concept. 

If  the  complex  field  on  a  line  is  Fourier  transformed,  one  can  identify  each 
Fourier  component  as  a  single  unique  plane  wave.  The  field  at  any  other  point 
can  then  be  calculated  by  adding  a  phase  shift  to  each  plane  wave  and  then 
summing  their  contribution  at  that  point.  In  an  attenuating  media,  this 
formulation  changes  slightly  because  the  phase  shift  includes  a  complex 
component  to  represent  the  attenuation. 

In  a  lossy  medium,  the  wavenumber  assumes  a  complex  value,  the 
imaginary  part  of  which  represents  the  effect  of  attenuation  in  the  medium. 
We  will  use  the  notation  kQ  for  the  wavenumber  in  a  lossless  medium  and  kc0 
when  the  medium  is  attenuating.  kco  is  defined  as: 

kco  =  K  +  J  fo  (2-1) 

If  the  field  energy  is  attenuated  in  the  medium  by  <;  db/cm  (the  unit  mostly 
used  in  microwave  imaging),  (0  would  be  linearly  related  to  <;  through  the 
following  relation, 

_  In  10 

£„  -  -  c  in  units  of  cm  .  (2.2) 
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u(x2ly)  =  /  A(ky,x2)  eJ  k»  y  dv 

27r  -CO 


(2.5) 


Also  it  must  satisfy  the  following  Helmholtz  equation 
V2u  +  kco2  u  =  0 


(2.6) 


Substituting  the  expression  in  Eq.  (2.5)  for  u  in  Eq.  (2.6),  we  get 


/ 


d2A(ky,x) 


dx‘ 


:=X2  -  ky2A(krx2)  4-  kc02  A(ky,x2) 


k’ y  dv  =  0  . 


(2.7) 


For  the  above  equation  to  be  satisfied,  the  term  inside  the  integral  has  to  be 
equal  to  zero. 


d2A(ky,x) 

dx2 


:=x,-(kco2  “  kv2)  -Mk  X2) 


(2.8) 


A  solution  for  A(ky,x2)  is  given  by 

A(ky,x2)  =  A(krx,)  ej  (2.9) 

This  equation  shows  the  exact  relationship  between  the  angular  spectra  of  the 
fields  measured  on  two  lines  separated  by  the  distance  (x2~x,).  This  is  exactly 
the  same  expression  for  the  propagation  of  the  angular  spectrum  in  a  lossless 
medium  [Gcodman68],  except  that  kc0  is  complex  in  this  case  and  thus  each 
Fourier  component  sees  both  a  phase  and  an  amplitude  change.  Substituting 
Eq.  (2.9)  in  Eq.  (2.5),  we  have 

u(x„.y)  =  —  7  A(ky,x.)  eJ  +  k,y|  dk  (2.10) 

*  2tt  _  } 


Now  recall  that  the  equation  for  a  unit-amplitude  plane  wave  propagating 
with  direction  cosines  (a.J)  is  simply 


i7l<>«  +  d  y  ) 
B(x,y)  =  o  x 


(2.11) 


where  X  is  the  wavelength.  Thus,  when  kro  is  real,  the  complex  exponential 
function  inside  the  integral  in  Eq.  (2.10)  may  be  regarded  as  a  plane  wave 
propagating  with  direction  cosines  (a.i).  where 


o 


(2.12a) 


and 
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(2.12b) 


When  kco  takes  complex  values,  a,  as  given  in  Eq.  (2.12a),  also  becomes 
complex.  Substituting  a  complex  value  of  a  in  Eq.  (2.11).  causes  each  of  the 
plane  wave  components  to  be  attenuated  in  the  x  direction.  Thus,  while  the 
direction  of  propagation  is  determined  by  the  direction  cosines 
(real  part  of  a  ,  J  ),  the  direction  of  attenuation  is  fixed  in  the  x  direction. 
Note  that  these  non-uniform  plane  waves  do  satisfy  the  homogeneous 
Helmholtz  equation  (  Eq.  (2.6)  ). 

After  the  above  reasoning,  the  following  question  may  arise.  Why  not 
consider  a  different  plane  wave  expansion  which  is  physically  more  realistic, 
namely  a  summation  of  plane  waves  for  which  the  directions  of  propagation 
and  attenuation  are  the  same?  For  such  an  expansion  ky  in  Eq.  (2.3)  has  to  be 
replaced  by  kcosin0,  where  0  is  the  angle  between  the  direction  of  propagation 
of  the  plane  wave  and  the  x  axis.  If  the  expansion  exists,  the  field  on  the  line 
x  =x,  could  be  expressed  as  (change  of  variable  kv  by  0 ) 

<X> 

u(x,.y)  =  —  /  A(^,x1)ejkoSin(<’)y  e'(»sin,#ly  kco  cos(0)  d0  .  (2.1-3) 
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The  last  exponential  term  inside  the  above  integral  will  blow  up  either  for  large 
negative  or  positive  values  of  y  depending  on  the  sign  of  6.  This  means  that 
the  field  on  the  line  x=x,  can  not  be  represented  as  a  summation  of 
attenuating  plane  waves. 


jc.v. 


2.2  Scattering  from  a  Single  Object 

In  general,  it  is  not  possible  to  calculate  a  closed  form  expression  for  the 
scattered  fields  from  arbitrarily  shaped  objects,  and  the  fields  must  be 
calculated  by  means  of  numerical  methods  such  as  finite  elements.  Usually,  to 
find  a  closed  form  solution  to  the  problem,  the  scattering  process  is  represented 
as  a  boundary  value  problem.  The  field  is  expressed  in  a  coordinate  system 
matching  the  boundaries  of  the  pioblern.  and  then  the  partial  differential 
equation  is  solved  by  applying  the  eigenfunction  expansions.  The  solution  to 
the  spattering  problem  is  much  simpler  and  easier  to  calculate  numerically 
when  the  object  is  considered  to  be  of  cylindrical  shape  rather  than  of  a  more 
complicated  shape. 

There  are  two  types  of  illumination  used  in  Diffraction  Tomograph) 
techniques:  a  plane  wave  or  a  point  source. 
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2.2.1  Plane  Wave  Illumination 

Assume  a  plane  wave  propagating  in  the  direction  x  is  incident  on  a 
cylinder  with  the  density  pe  and  the  compressibility  Ke  centered  at  the  origin, 
embedded  in  a  medium  with  a  density  of  p  and  a  compressibility  of  /c  (see  Fig. 
2.1).  The  complex  envelope  of  the  incident  field,  defined  as  the  field  present 
when  the  cylinder  is  removed  from  the  medium,  is  expressed  as 

Uj(x.y)  =  U0  ej  k°°x  (2.14a) 

or  in  polar  coordinates  as 

u,(T)  =  U„  ej  k~ r  cos(*)  (2.14b) 

where  kf0  is  the  wavenumber  and  (r,0)  are  the  coordinates  of  the  point  (x,y) 
defined  in  cylindrical  coordinates.  The  cylindrical  boundaries  of  the  object, 
consisting  of  a  single  cylinder,  suggest  the  eigenfunctions  to  be  of  cylindrical 
type.  Therefore,  the  incident  plane  wave  should  be  expressed  in  terms  of  these 
cylindrical  eigenfunctions  as  [Morse68,  \Veeks64] 

'ii(T)  =  U0  V  £n  j"  cos(no)  Jn(kf0r)  ,  (2.15) 

n  =0 

where  £n  is  1  for  n=0  and  2  otherwise  and  Jn(.)  is  the  nth  order  Bessel  function. 
When  the  cylinder  with  its  axis  at  r=0  is  present,  the  field  is  no  longer  given 
by  either  of  Eqs.  (2.14a)  or  (2.14b).  Instead  there  is  present,  in  addition  to  the 
incident  plane  wave,  a  scattered  outgoing  wave.  In  general  the  form  of  this 
scattered  wave  is  given  by  the  series  [Morse68] 


"5m  = 


v  An  cos(no)  Jn(kcr) 

n  =0 


V  Bn  cos(nd>)  j  .Jn(kror)  +  j  Nn(kP0r) 

n  =0 


r  <  a 


r  >  a 


(2.16) 


where  k,,  is  the  wavenumber  inside  the  object  and  Nn(.)  is  the  nth  order 
Neumann  function.  The  combination  I+jN.  used  for  r  >  a,  is  to  ensure  that 
all  of  the  scattered  energy  travels  away  from  the  cylinder,  and  J,  used  for 
r  <  a.  is  to  ensure  that  the  scattered  energy  travels  toward  the  interior  of  the 
cylinder.  The  boundary  conditions  to  be  imposed  on  the  problem  are  the 
continuity  of  the  pro-sure  and  the  axial  velocity  across  the  object-medium 
interface  for  the  case  of  ultrasound  energy. 
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Figure  2  1 

The  scattered  field  from  a  cylinder  of  radius  a.  illuminated  by 
wave,  is  measured  at  the  point  [r ,<j>). 


a  plane 
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In  a  visco-elastic  rr.edium,  the  complex  envelope  of  the  radial  velocitv,  V 
and  the  pressure,  P,  are  related  by  the  following  relationship 
1  <9P 

J  u  Vr  = - —  ,  (2.1/ » 

p  dr 

where  is  the  angular  frequency  of  the  incident  field.  We  take  us,  defined  as 
the  scattered  field  in  Eq.  (2.16),  to  be  the  scattered  pressure  field.  By  imposing 
the  boundary  conditions  (continuity  of  Vr  and  P)  at  the  surface  of  the  cylinder, 
one  can  solve  for  the  coefficients  An  and  Bn,  defined  in  Eq.  (2.16).  Since  we  are 
only  interested  in  the  scattered  fields  outside  the  object,  we  only  bring  the  part 
of  solution  for  r  >  a,  given  by  [Morse68j 


us01  -  S  4  jn  Dn  cos(n?)  Hnl»(kj7 


(2.18) 


where 


Dn  = 


Jn(^ea)  J' n(^coa)  a  ^  n(^ea)  ^nt^co3-) 

aJ'n(kea)  Hn“'(keoa)  -  Jn(kea)  H„*l)'(kcoa) 


(2.10) 


o-  =  —  ,  (2.20) 

K  pe 

and  £n  is  1  for  n=0  and  2  otherwise,  Jn(.)  is  an  nth  order  Bessel  function,  H^1 '( . ) 
is  a  Hankel  function  of  the  nlh  order,  kco  and  ke  are,  respectively,  the 
wavenumbers  in  the  medium  and  the  object,  a  is  the  radius  of  the  cylinder,  0  is 
the  angle  between  the  direction  of  the  incident  plane  wave  and  the  line  joining 
the  center  of  the  cylinder  and  the  observation  point,  and,  finally,  p  and  k  are 
the  density  and  the  compressibility  constants  of  the  medium. 

When  the  object  is  lossy,  ke  assumes  a  complex  value  and  only  the 
coefficients  Dn  in  Eq.  (2.18)  will  change.  However,  when  the  medium  is 
attenuating  (complex  kco),  the  Hankel  function  in  Eq.  (2.18)  is  affected  as  well. 
Because  of  the  attenuation  the  Hankel  function  dies  out  exponentially  as  a 
function  of  |t|  .  This  is  clear  by  considering  the  asymptotic  form  of  this 
function, 

i  j  (  z  -  n- ~  1  ') 

'(z)  -  -e  -  (2.21) 

z 

This  behavior  agrees  with  <>nr  physical  intuition  of  the  problem,  since  th“ 
energy  from  a  source  should  decay  rapidly  as  a  function  of  distance. 
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2.2.2  Point  Source  Illumination 

In  some  reconstruction  techniques,  instead  of  illuminating  the  object  with 
a  plane  wave  a  small  sized  transducer  is  used  as  the  transmitter.  We  will  talk 
about  these  reconstruction  techniques  in  later  chapters.  A  transducer  with  a 
very  small  sized  aperture  can  be  approximated  by  a  point  source. 

Assume  a  point  source  located  at  70  is  illuminating  a  cylinder  with  the 
density  pe  and  the  compressibility  Ke  embedded  in  a  medium  with  the  density  p 
and  the  compressibility  k  (see  Fig.  2.2).  The  complex  envelope  of  the  incident 
field  at  7  is  expressed  by 

u,(7)  =  UoHo01(kco|7-7o|  )  .  (2.22) 

where  kco  is  the  wavenumber.  The  cylindrical  boundaries  of  the  object, 
consisting  of  a  single  cylinder,  suggest  the  eigenfunctions  to  be  of  cylindrical 
type.  Using  the  expansion  of  the  Hank'd  function  in  cylindrical  coordinates 
[Morse53j,  we  express  the  incident  wave  in  the  following  form 


Uj(7)  -  U0  V)  ‘n  cos  n(0-<po) 

n  =  0  1 


Jn(kcor)  Hn(1)(kcoro)  r  <  T0 


Jn(kcor0)  Hn0)(kcor) 


r  >  r0- 


(2.23) 


When  the  cylinder  with  its  axis  at  the  origin  is  present  the  total  field  would  be 
the  incident  w'ave  plus  the  scattered  wave.  In  general  the  form  of  this 
scattered  wave  is  given  by  the  series  in  Eq.  (2.16).  Again,  the  boundary 
conditions  imposed  on  the  problem  are  the  continuity  of  the  pressure  and  the 
axial  velocity  across  the  object-medium  interface. 

By  imposing  the  boundary  conditions  at  the  surface  of  the  cylinder,  we 
solve  for  the  coefficients  An  and  Bn.  Since  we  are  only  interested  in  the 
scattered  field  outside  the  object,  we  only  bring  the  result  for  r  >  a. 


>i5m  -  1]  Jn  Dn  cos  n(0-d>o)J  Hn,n(k,0|  r|  ) 

n  =0 


where 


D„  = 


Hn<M  rol  1  MM)  I'n(kr„a)  -  J'n(M)  In(kroa) 
«J'n(M)  Hn(,,(ka)  -  Jn(M)  Hn"i'(ka) 


where  the  parameters  are  the  same  as  those  defined  in  section  2.2.1. 
Considering  Eqs.  (2.  IS).  (2.19),  (2.21)  and  (2.23),  we  notice  that  the  coefficients 
I)n  are  different  by  a  factor  of  IIn(k.„|7„|  )  and  also  the  argument  of  the  cosine 
fiinct  ion  m  I  ■  |  (2  21)  includes  the  posit  jonal  angle  of  the  point  source  at  7(1. 
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If  the  medium  is  attenuating,  the  same  effects  as  we  discussed  in  section 
2.2.1  would  be  observed.  As  the  point  source  at  70  is  moved  away  from  the 
cylinder  the  term  Hn(kco|70|  )  in  Eq.  (2.25)  will  decrease  and  the  coefficients  Dn 
would  get  smaller  and  so  would  the  magnitude  of  the  scattered  field,  as  is 
expected. 

2.2.3  Numerical  Calculation  of  Scattered  Fields 

We  have  based  our  computer  simulation  programs  on  the  expressions 
derived  in  Eqs.  (2.18)  and  (2.24).  The  field  is  calculated  on  a  straight  receiver 
line  on  the  opposite  side  of  the  object  relative  to  the  direction  of  illumination. 
All  Diffraction  Tomographic  techniques  in  transmission  mode  require  such  a 
geometry  for  the  receiver  positions.  Since  the  amount  of  data  required  for 
reconstructing  the  cross  section  of  an  object  is  relatively  large,  it  is  beneficial  to 
generate  efficient  code  to  speed  up  the  computation.  We  have  used  the  array 
processor  FPS  AP120B  to  carry  out  the  vectorized  operations.  Most  of  the 
computational  effort  in  our  problem  is  spent  on  the  calculation  of  the  Bessel 
functions,  Jn,  and  the  Neumann  functions,  Nn.  A  Neumann  function  is  the 
imaginary  part  of  a  Hankel  function  of  the  first  kind,  Hn^. 

When  both  the  object  and  the  medium  are  lossless,  we  calculate  the 
functions  J0,  J1(  N0  and  Nj  by  means  of  non-vectorized  codes  available  in  the 
UNIX  mathematical  library.  Higher  order  terms  can  be  calculated  through  the 
following  recursive  formulas  [Bowman58] 

Jn  +  I(x)  =  J„(x)  ,  (2.26) 

and 

N„  +  ,(x)  =  N„(x)-N„.,(x)  .  (2.27) 

By  forward  recursion,  with  correct  values  of  J0(x)  and  J,(x).  one  fails  to 
produce  satisfactory  values  for  Jn(x)  much  beyond  n  =  x.  This  is  due  to  the 
large  factor  multiplied  by  Jn(x)  in  Eq.  (2.26),  which  causes  any  initial 
uncertainties  and  rounding  errors  to  propagate  in  multiples  of  that  factor. 
Knowing  that  Jn(x)  goes  to  zero  as  n  goes  to  infinity,  one  can  reformulate  the 
problem.  Starting  with  Jn  +  m(x)  =  0  and  the  arbitrary  Jn  +  m_,(x)  =  1,  for  some 
sufficiently  large  m,  the  solution  for  .J n ( x )  is  sought  by  backward  recurrence 
[Fox 68],  With  a  "check"  value  at  some  n,  say  n=0,  the  '‘trial"  solution  can  be 
scaled  to  satisfy  this  condition.  It  should  be  noted  that  recursive  calculation  of 
the  Neumann  function  does  not  suffer  from  this  complication.  VV  e  have 


developed  vectorized  programs  for  these  recursive  calculations  (forward 
recursion  as  well  as  backward  recursion),  and  in  fact,  calculation  of  all  the 
required  higher  order  terms  takes  almost  no  time  in  comparison  with  that  of 
zeroth  and  first  order  terms  which  are  calculated  by  the  non-vectorized  code. 

The  derivative  of  the  Bessel  functions  and  the  Neumann  functions  are 
calculated  through  the  following  relations  [Bowman58j: 

J'„W  =  {(•>„-.(*)  -•>„*.(*)]  •  (2  28) 

and 

N'.M  =  7  (  Nn_,(x)  -  N„  +  ,(x)  )  .  (2  29) 

When  the  object  is  lossy  (ke  is  complex),  by  considering  equations  (2.18) 
and  (2.24),  one  realizes  that  only  the  coefficients  Dn  will  be  different  from  the 
case  of  a  lossless  object.  This  means  that  for  all  possible  positions  of  the 
receiver,  the  coefficients  Dn  are  fixed  and  need  to  be  calculated  only  once  for 
each  cylinder.  Therefore  the  computation  time  for  the  calculation  of  scattered 
fields  from  lossy  objects  would  be  on  the  same  order  as  that  of  the  lossless 
objects.  According  to  Eqs.  (2.19)  and  (2.25),  calculation  of  the  coefficients  Dn, 
requires  a  routine  providing  Bessel  functions  of  complex  arguments,  so  we  made 
use  of  the  “mmbzjn”  routine  in  the  well  known  “IMSL”  library. 

When  both  the  object  and  the  medium  are  lossy,  the  computational 
difficulties  increase  severely.  This  time  kco,  the  wavenumber  in  the  medium 
becomes  a  complex  number.  This  affects  the  calculation  of  coefficients  Dn, 
because  not  only  Bessel  functions  with  complex  arguments  are  required,  as  in 
the  case  of  lossy  object,  but  also  Neumann  functions  with  complex  arguments 
are  needed  as  well.  Unfortunately  we  were  not  able  to  find  a  routine  that 
calculates  the  Bessel  functions  of  the  second  kind  in  any  of  the  software 
libraries  at  Purdue.  However,  after  a  long  search  for  computational  procedures 
for  calculation  of  Neumann  functions  with  complex  arguments,  we  found 
several  useful  references  [Goldstein59,  Marion61,  StegunoT].  Finally,  we  were 
led  to  the  following  procedure  for  obtaining  Bessel  functions  of  the  second  kind 
with  complex  arguments. 

To  compute  N0(z),  the  zeroth  order  Neumann  functions  with  complex 
arguments,  we  use  the  following  series 
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N0(z)  =  f  [  C  +  log(f)  J  J0(z)  +  2  V  (~1  )n  +  *  ±  J2n(z)  .  (2.30) 

l  2  n=l  n 


In  this  equation  the  constant  C  is  the  Euler's  constant.  N,(z)  can  be  found 
using  the  relation 

J,(z)  N0(z)  -  J0(z)  N,(z)  =  ~  .  (2.31) 

The  functions  Nn(z),  n=2,3,..  may  now  be  obtained  from  the  usual  recurrence 
relation. 


Np  +  i(z)  =  2  p  z  Np(z)  -  Np_,(z)  .  (2.32) 

Although  the  calculation  of  N0  is  very  expensive,  higher  order  terms  are  simply 
calculated  from  Eqs.  (2.31)  and  (2.32).  Thus,  in  principle,  the  calculation  of 
scattered  fields  for  attenuating  objects  and  media  is  not  much  more  expensive 
than  that  for  the  case  of  lossless  material,  when  non-vectorized  code  is  used  in 
both  cases. 


As  mentioned,  in  the  case  of  a  lossless  medium,  our  computer  program  is 
developed  so  that  the  zeroth  and  the  first  order  Hankel  functions  are  calculated 
using  non-vectorized  routines  and  then  fed  to  the  FPS-AP120B  array  processor 
for  vectorized  calculation  of  the  higher  order  terms.  An  attempt  to  implement 
a  similar  program  for  Bessel  functions  with  complex  arguments  failed  due  to 
the  complex  nature  of  the  algorithm.  With  more  work,  the  calculation  of 
Neumann  functions  could  be  vectorized  but  since  it  can  only  speed  up  the 
calculation  by  a  factor  of  two  (the  calculation  of  Jn(z)  could  not  be  vectorized) 
we  do  not  think  it  is  justifiable. 

The  computer  program  developed  was  used  to  generate  the  scattered  fields 
from  a  cylinder  with  diameter  of  6  wavelengths  illuminated  by  a  5  MHz  plane 
wave.  In  Figs.  2.3  and  2.4,  we  have  shown  the  magnitude  of  the  scattered  field 
from  a  lossless  cylinder  with  the  solid  line  and  that  for  a  lossy  cylinder  with  the 
dashed  line.  The  object's  refractive  index  is  of  1.1  and  1.5,  respectively.  In 
both  cases  the  cylinder  is  embedded  in  a  lossless  material  and  the  receiver  line 
is  10  wavelengths  from  the  center  of  the  cylinder.  Considering  these  two 
figures,  one  observes  that  while  the  field  magnitude  decreases  in  the  presence  of 
the  attenuation  in  the  object,  the  overall  structure  of  the  scattered  field  does 
not  change  appreciably. 

In  Fig.  2.5,  we  have  shown  the  efToct  of  an  attenuating  medium  on  the 
scattered  fields  from  an  object.  The  solid  line  in  Fig.  2.5  represents  the 
magnitude  of  the  scattered  field  from  a  lossless  cylinder  with  diameter  of  G 
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Figure  2.3 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.1  (embedded  in  a  lossless  medium)  illuminated  by  a  plane 
wave,  measured  on  the  receiver  line  at  a  distance  of  10  X  from  the 
cylinder.  The  solid  line  corresponds  to  the  case  of  a  lossless  cylinder,  while 
the  dashed  line  represents  the  field  when  the  attenuation  constant  of  the 
cylinder  is  3.5  db/cm. 


•  V 

0-J 

\y 

*4  „  - 


position  on  the  receiver  line  <  *  X  > 


Figure  2.4 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.5  (embedded  in  a  lossless  medium)  illuminated  by  a  plane 
wave,  measured  on  the  receiver  line  at  a  distance  of  10  X  from  the 
cylinder.  The  solid  line  corresponds  to  the  case  of  a  lossless  cylinder,  while 
the  dashed  line  represents  the  field  when  the  attenuation  constant  of  the 
cylinder  is  3.5  db/cm. 
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Figure  2.5 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.1  illuminated  by  a  plane  wave,  measured  on  the  receiver  line 
at  a  distance  of  10  X  from  the  cylinder.  The  solid  line  corresponds  to  the 
case  of  a  lossless  cylinder  and  medium,  while  the  dashed  line  represents 
the  field  when  the  attenuation  constants  of  the  cylinder  and  the  medium 
are  both  2.5  db/cm. 


wavelengths  and  refractive  index  ef  1.1  embedded  in  a  lossless  medium  while 
the  dashed  line  corresponds  to  the  case  where  the  attenuation  constant  of  the 
medium  and  the  cylinder  are  both  2.5  db/cm.  In  each  case  an  attenuating 
plane  wave  was  incident  on  the  cylinder  and  the  field  was  measured  aiong  the 
receiver  line,  10  wavelengths  away  from  the  cylinder.  Fig.  2.6  reflects  the  same 
information  for  a  larger  refractive  index  of  1.5.  These  figures  show  that  not 
oulj  the  attenuation  causes  a  iarge  loss  of  energy  but  it  also  causes  the  field  on 
the  receiver  line  to  be  multiplied  by  a  smooth,  bell-shaped  window.  This 
window  is  due  to  the  direct  relationship  between  the  attenuation  of  the  field 
and  the  range  of  propagation. 

Figs  2.7-2.10  show  the  same  information  as  in  Figs  2.3-2. 6,  while  the 
incident  plane  wave  is  replaced  with  a  point  source  10  wavelengths  away  from 
the  cylinder. 

2.3  Scattering  from  a  Multi-component  Object 

One  can,  in  principle,  obtain  the  exact  solution  to  the  wave  equation  for  a 
multi-component  object  provided  one  is  able  to  solve  the  boundary  value 
problem  for  the  entire  object.  In  practice  it  is  not  possible  to  do  so  even  for 
two-  or  three-component  objects,  and  one  must  take  recourse  to  computational 
procedures. 

For  calculating  the  scattered  fields,  a  major  source  of  difficulty  with 
multi-component  objects  is  dealing  with  the  interaction  between  the  various 
components.  Depending  on  the  interaction  between  the  components,  the  total 
scattered  field  may  or  may  not  bear  any  resemblance  to  the  simple  sum  of  the 
scattered  fields  for  each  of  the  components,  assuming  the  others  to  be  absent. 
In  this  section  we  have  presented  a  new  computational  procedure  for 
calculating  the  inter-component  interaction.  With  the  computer  programs 
developed  to  date,  we  now’  are  able  to  generate  the  scattered  fields  that  are  not 
limited  by  the  first  order  assumption.  Although  the  results  shown  will  only 
include  the  second  order  fields  for  a  multi-component  object,  the  computational 
procedure  can  easily  be  generalized  for  higher  order  scattering  effects. 

In  what  follows,  we  have  first  briefly  reviewed  the  basic  theory  of  multiple 
scattering  (Section  2.3.1).  This  is  followed  by  t  ho  discussion  of  the 
computational  procedure  for  the  calculation  of  the  multiple  scattering  effects 
(Section  2.3.2).  Results  on  two  component  objects  are  shown  in  Section  2  3  3 
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Figure  2.8 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.5  (embedded  in  a  lossless  medium)  illuminated  by  a  point 
source  10  X  away  from  the  cylinder,  measured  on  the  receiver  line  at  a  dis¬ 
tance  of  10  X  from  the  cylinder.  The  solid  line  corresponds  tc  the  case  of 
a  lossless  cylinder,  while  the  dashed  line  represents  the  field  when  the  at¬ 
tenuation  constant  of  the  cylinder  is  3.5  db/cm. 
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Figure  2.9 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.1  (embedded  in  a  lossy  medium)  illuminated  by  a  point  source 
10  X  away  from  the  cylinder,  measured  on  the  receiver  line  at  a  distance  of 
10  X  from  the  cylinder.  The  solid  line  corresponds  to  the  case  of  a  lossless 
cylinder,  while  the  dashed  line  represents  the  field  when  the  attenuation 
constants  of  the  cylinder  and  the  medium  are  both  2.5  db/cm. 
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Figure  2.10 

Magnitude  of  the  scattered  field  from  a  cylinder  of  radius  3  X  and  refrac¬ 
tive  index  1.5  (embedded  in  a  lossy  medium)  illuminated  by  a  point  source 
10  X  away  from  the  cylinder,  measured  on  the  receiver  line  at  a  distance  of 
10  X  from  the  cylinder.  The  solid  line  corresponds  to  the  case  of  a  lossless 
cylinder,  while  the  dashed  line  represents  the  field  when  the  attenuation 
constants  of  the  cylinder  and  the  medium  are  both  2.5  db/cm. 
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2.3.1  Review  of  the  Basic  Formulation  of  Multiple  Scattering 

We  have  based  our  algorithms  on  basic  theory  in  which  the  scattered  field 
is  expressed  as  an  infinite  summation  of  different  order  terms  [Ishimaru78j. 
First  order  fields  are  obtained  by  considering  the  interaction  of  the  original 
incident  field  with  each  component,  assuming  the  others  to  be  absent.  First 
order  fields  caused  by  one  component  when  incident  on  the  other  components 
generate  the  second  order  fields,  and  so  on. 

The  b  asic  elements  of  multiple  scattering  theory  that  we  will  use  in  our 
computational  modeling  will  now  be  explained  in  detail  with  the  help  of  a  two 
component  object  (The  approach  is  easily  generalized  to  more  than  two 
components).  Fig.  2.11  shows  a  two  component  phantom  being  illuminated  by 
a  field  denoted  by  U;(7)  (In  the  absence  of  the  phantom  the  field  everywhere  will 
be  Ujfr}).  We  will  use  u(7a)  to  denote  the  actual  field  at  a  point  frj  as  shown  in 
the  figure.  u(Ta)  is  equal  to 

u(TJ  =  u,fraj  +  0,fr a)  +  p,fTa)  .  (2.33) 

where  dj(7a)  and  <^>2(7a)  are,  respectively,  the  scattered  fields  at  7a  caused  by  the 
phantom  components  1  and  2. 

Let  uinc(T;)  be  the  incident  field  at  the  center  of  the  ith  component,  and  let 
°s(7a7j)  be  an  operator  function  such  that  when  it  is  applied  to  the  fields 
incident  ^n  the  scattering  component  at  it  generates  the  scattered  field  at 


the  observation  point  7a.  In  terms  of  o3’s,  0,  and  <j>2  are  given  by 

<*i(Ta)  =  os(rvr,)uincffi)  i  =  L2  (2.31) 

Substituting  (2.31)  in  (2.33).  we  get 

uCFJ  =  u.ffj  +  o?ffa|7,)uinc(r,)  +  o,(ra|r2)uinc(To)  (2.35) 

The  field  incident  at  the  site  of  each  scatterer  may  be  expressed  as 

uinc(ri)  =  ui(F|)  +  o#rr,|72)uinc(7.)  .  (2.36a) 

uin,.(T2)  =  U|fr2)  +  o  (T>jT,)uin.(T|)  (2.36b) 

Substituting  (2.36a,b)  in  (2.35).  we  get 

h(T  )  ~  u  fT\)  <>  (T.J  7|  )u  rf, )  <>  (7aj  72)u;(7..) 

+  "K  T,  K(7,  i  72)ui:,..(r2)  ^  ".rra;  72).'  (72j  T,  )u (T*, )  (2.371 


If  we  again  substitute  (2  36a. b)  in  (2.37).  ue  have  the  following  expression 


Ou^ 


Incident 

field 


component  I 


Figure  2.11 

A  two  component  phantom  illuminated  with  incident  field  w  (F) 
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u(TJ  =  u,(Tj  +  osfTa|r,)ui(r1)  +  osf7a|  FoKfTo) 

+  os0ra|'ri)Osfri|f2)u1(T2)  +  os(7a  j  72  )os(72  j  7, )  u  ,(7, ) 

+  °s(TaJ  Tj)os  frj  j  72)os(72  j  T,  )iiinc(T, ) 

+  o8(?;|T2)os(r2|T|)o9(Fi|T2)uinc(F3)  (2.38) 

A  more  compact  way  to  express  the  above  equation  is 

U(TJ  =  u,(fa)  +  ufir5t  -orderf^a^  ^second-order(Ta) 

+  Uhigher-order(Ta)  (2.39) 

The  quantities  ufirst_order,  usecond_order  and  uhigher_order  represent  the  first  order, 
second  order,  and  higher  order  contributions  at  the  observation  point  Ta.  These 
are  given  by 

Ufirst-orderfa)  =  °,(T» |  71 ) « i(T, )  +  Osf7a|  72)U;(7o)  (2.40a) 

^second-orderC^a)  =  To)os(72 1 T,  )uifF1 ) 

+  Osteal  7,)og(7,  I  T2)u,(T2)  (2.40b) 

Uhigher-ordertTa)  —  ^sC^a  |  r  1  |  |  r2  (  rl  ) 

+  Os(Ta|r2)os(r2|T,)os(T,|r2)uinc(r2)  (2.40c) 

In  Fig.  2.12,  we  have  shown  the  first  and  the  second  order  scattering  processes 
for  a  two  component  object.  For  the  first  order  scattering,  each  component 
interacts  independently  with  the  incident  illumination,  being  oblivious  of  the 
existence  of  the  other.  To  compute  the  second  order  scattering  terms,  each 
component  interacts  with  the  fields  sent  in  its  direction  by  the  other 
component,  and  so  on. 

To  compute  the  second  term  in  (2.39),  which  is  given  by  equation  (2.40a). 
we  require  o.(T3|Tj)Uj(7j),  where  7,  is  the  location  of  the  scatterer  and  7a  the 
observation  point.  When  the  components  are  cylinders  as  in  Fig.  2.12.  and 
when  the  illumination  field  is  a  plane  wave,  the  operator  is  given  by  the 
expression  in  Fq  '2 .1^).  When  the  source  of  illumination  is  a  point  source,  one 
should  consider  t  expression  in  Fap  (2  21)  instead. 


Receivers 

line 


Receiver 

line 


Figure  2.12 

Th  is  figure  depicts  the  first-  and  the  second-order  scattered  fields  for  a  two 
component  object. 
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2.3.2  A  New  Procedure  for  Computing  Multiple  Scattering 

In  this  section  we  will  use  the  notation  k0  for  the  wavenumber  in  the 
medium,  and  assume  that  there  is  no  attenuation  present  in  the  medium.  The 
reason  will  be  explained  in  a  later  section  (2.3.4).  We  can  directly  use  Eq.  (2.18) 
or  (2.24)  to  compute  the  first  order  scattering  term  in  (2.40a);  it  is  the 
computation  of  the  second  and  higher  order  terms  that  poses  difficulty.  The 
computational  procedure  presented  in  this  section  addresses  the  calculation  of 
the  second  order  term  as  given  by  (2.40b).  Generalization  of  this  procedure  for 
computing  the  higher  order  scattering  effects  is  straightforward. 

By  using  (2.18)  or  (2.24),  w'e  first  calculate  the  scattered  fields  as  generated 
by  the  interaction  of  the  incident  field  with  the  first  component  on  the  line  AV 
as  shown  in  Fig.  2.13.  For  this  purpose,  we  choose  a  new  set  of  Cartesian 
coordinates  denoted  by  (£ ,rj )  such  that  the  r)  axis  is  along  the  vector  from  7t  to 
72.  The  origin  of  the  new  coordinate  system  is  at  the  midpoint  between  7,  and 
72;  the  position  vector  for  the  origin  is  therefore  given  by 


7,  +  72 


(2.41) 


We  will  use  A(a)  to  denote  the  Fourier  transform  of  the  scattered  field  on  this 
line.  It  is  given  by 


A(o)  -  /  os(£,f/=0|7,)  Uj(7,)  e“JO*  d£  (2.42) 
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In  terms  of  the  Fourier  transform  A(rv),  the  field  at  any  point  in  the  half  space 
r]>0  can  be  expressed  as 

os(7j7,)  11,(7, )  =  f  A(o)  ejRK'  {f  7n)  do  (2.43) 

2  K 

where  K'  —  (o,  \/k0~- a*)  and  R  is  the  rotational  transformation  applied  to  the 
rj  coordinates  to  position  it  along  the  x—  y  coordinates  (see  Fig.  2.14).  N  oto 
that  o  is  the  £  component  of  f\  ,  and  ^/k,,"-o_>  the  component. 

When  the  magnitude  of  o  exceeds  the  wavenumber  k  .  the  component  of 
hd  in  the  r/  direction  becomes  imaginary.  'The  corresponding  wave  components, 
called  evanescent  waves,  will  be  attenuated  as  they  travel  in  rj  direction:  the 
larger  the  value  of  rv.  the  greater  the  attenuation  Beyond  a  distance  of  about 
10  wavelengths  from  the  source,  the  evanescent  waves  are  usually  unimportant. 


/.v'/.V. 


/  u 


cylinder  2 


cylinder  1 


Figure  2.14 

A  new  set  of  Cartesian  coordinates  (£,»7)  is  chosen  such  that  the  r;-axis  is 
along  the  vector  from  7X  to  To  and  the  origin  is  at  the  midpoint  between 
T,  and  To.  ^ '  *s  wavenumber  vector  in  the  new  coordinates  and 
(7-7 0 )  is  the  vector  position  of  the  point  7  in  the  new  coordinates. 
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For  digital  implementation,  we  rewrite  equations  (2.42)  and  (2.43)  in  the 
following  discrete  form 


T-' 


Am  =  A£  V)  os(nA£,  t]~ 0  |  T,)  u,(7,)  e  N 
_  N 


JTT  m  n 


(2.44) 


?i)  UifF,)  - 
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A  e1 


jRK' m  .  (?-?„] 


f" 


(2.45) 


where  A£  is  the  sampling  interval  in  the  £  direction,  and  N  is  the  total  number 
of  sampling  points.  The  variable  R' m  is  given  by 


,  £7mii  /.  o  ,  2 trm  ,  *  .  . 

V  ~  NA^  ’  V  k  NA£  "'46 

The  operation  in  (2.44)  can  be  implemented  as  a  Fast  Fourier  Transform 
(FFT).  When  the  argument  of  the  square-root  becomes  negative,  the 
corresponding  wave  becomes  evanescent.  This  evanescent  wave  propagates  in 


the  £  direction  with  a  spatial  frequency  of 


-,  while  it  undergoes 


attenuation  by  the  factor  exp  ~r/- 


•)2~k02  in  the  r/  direction.  The 


higher  the  spatial  frequency  in  the  £  direction  is,  the  greater  the  attenuation  in 
the  i /  direction  will  be. 

The  sampling  interval  is  determined  by  the  choice  of  the  extent  to  which 
on.-  wishes  to  include  the  evanescent  waves  in  the  measurement.  The  parameter 
7  will  be  a  measure  of  this  inclusion.  Its  value  will  range  from  0  to  1;  when  it  is 
equal  to  0  the  sampling  interval  will  be  such  as  to  exclude  the  evanescent 
waves.  When  it  is  equal  to  1,  all  evanescent  waves  will  be  included.  We  obtain 
A£  by  using  the  following  equality 


,  N  , 

1 

-•'71 

L  2 

[  NAC 

—  K  . 

e  1  k  '  =  7  12.17) 

The  left  hand  side  is  the  attenuation  suffered  by  the  evanescent  wave, 
characterized  by  the  highest  spatial  frequency  along  the  £  direction.  A£  can  be 
calculated  from  (2.47)  to  be 


s  s 


For  7  equal  to  1,  this  expression  yields  a  value  of  X/2  for  the  sampling  interval, 
a  result  that  is  intuitively  justifiable. 


It  is  less  possible  to  be  analytic  about  the  choice  of  X,  the  number  of  data 
points  on  the  £  axis.  The  interval  between  the  samples  in  the  angular 
spectrum  is  determined  by  the  total  number  of  sampling  points  N.  This 
sampling  interval,  and  therefore  the  value  of  X,  should  depend  on  how  fiat  or 
sharply  peaked  the  angular  spectrum  representation  is.  Ideally  the  best  choice 
for  X  is  that  for  which  no  undersampling  occurs  in  the  discrete  representation 
of  the  angular  spectrum. 


The  second  order  scattered  fields  on  the  receiver  line  may  then  be 
determined  by  applying  the  formula  in  (2.18)  and  (2.2-1)  to  each  plane  wave 
component.  Therefore, 


V  \  --iRK'm-ffj-M 
XA£  “'N  '  m 


x|expr.  in  eq.(2.18)  or  (2.21)  with  r  replaced  by  Toj  (2.49) 

Third  order  scattering  terms  can  be  calculated  in  the  same  way;  this  time  the 
second  order  fields  are  first  calculated  on  the  line  AA  instead  of  the  receiver 
line;  and  then  decomposed  into  plane  waves.  The  only  harrier  to  the 
calculation  of  higher  order  terms  is  the  computer  time  involved. 

Although  the  rate  of  convergences  of  the  series  in  Kqs.  (2.18)  and  (2.24) 
are  dependent  on  the  change  in  refractive  index,  it  is  much  more  sensitive  to 
the  value  of  real(kro)a,  the  higher  the  value  of  real(kro)a  is,  the  slower  the 
convergence.  In  all  the  computer  simulation  results  shown  in  this  paper,  the 
number  of  terms  retained  on  the  right  hand  side  in  Kqs.  (2.18)  and  (2.24)  was 
100.  this  number  being  more  than  adequate  to  render  the  truncation  errors 
negligible  in  all  cases  of  interest  to  us.  Although  we  could  have  used  fewer  than 
100  terms  without  sacrificing  accuracy  for  some  values  of  real(kro)a  and  other 
parameters,  due  to  the  programming  requirements  of  the  Floating  Point  Array 
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Processor  (AP120B)  that  was  used  to  generate  the  results,  we  preferred  using  a 
fixed  number  of  terms. 

If  we  only  compute  the  second  order  contributions  to  the  scattered  fields, 
the  computing  effort  required  is  proportional  to  0(KW12)  ,  where  K  is  the 
number  of  terms  retained  on  the  right  hand  side  in  Eq.  (2.18)  or  (2.24),  N  is  the 
number  of  plane  wave  components  in  the  expansion  of  the  fields  on  line  AA’ , 
and  M  the  number  of  object  components.  For  a  two-component  object,  the 
CPU  time  was  less  than  1  minute.  For  computing  higher  order  scattering 
contributions,  the  computation  time  is  proportional  to  0(KNMq)  ,  where  q  is 
the  maximum  order  of  scattering  included  in  the  simulation. 

2.3.3  Numerical  Results  on  the  Calculation  of  Multiple  Scattering 

Although  the  results  shown  will  only  include  the  second  order  fields  for  a 
multi-component  object,  the  computational  procedure  can  easily  be  generalized 
for  higher  order  scattering  effects. 

We  will  now  show  the  results  obtained  with  the  computer  implementation 
of  the  procedure  discussed  in  the  preceding  section.  In  this  section,  all  our 
results  are  on  two  component  objects,  and  only  the  scattering  contributions  up 
to  the  second  order  have  been  included.  When  the  ratio  of  back-scatter  to 
forward-scatter  in  a  two  component  object  i>  small,  which  happens  when  the 
components  are  large  compared  to  a  wavelength,  the  scattering  contributions  of 
orders  higher  than  two  become  negligible.  This  condition  is  satisfied  in  a 
majority  of  the  simulation  results  shown;  while  for  others,  computational 
expediency  was  our  only  reason  for  not  going  beyond  the  second  order  scatter. 
W  e  should  like  to  add  that  even  when  the  object  components  are  large,  it  may 
still  be  necessary  to  include  the  higher  order  scattering  terms  if  there  are  more 
than  two  components  involved. 

The  results  shown  in  Figs.  2.15  through  2.17  and  2.20  through  2.26  are  for 
different  situations  with  regard  to  whether  the  components  do  or  do  not  block 
each  other,  their  sizes,  and  the  deviations  in  their  refractive  index  values  from 
the  background.  Our  intention  is  not  to  inundate  the  reader  with  a  large 
number  of  plots.  The  plots  were  selected  either  bacause  they  illustrate  how 
rapidly  the  higher  order  scattering  contributions  become  important  as  the 
controlling  parameter  (refractive  index  deviation,  diameter,  etc.)  was  changed 
incrementally  over  one  or  two  steps,  or  because  they  represented  different 
situations  such  as  the  components  blocking  or  not  blocking  each  other,  etc. 
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There  are  four  plots  for  each  case  considered  in  Figs.  2.15  through  2.17 
and  Figs.  2.20  through  2.26.  Plots  labeled  (a)  show  the  scattering  geometry, 
such  as  the  position  of  the  cylinders  relative  to  the  receiver  line  (the  vertical 
line  on  the  right)  and  their  sizes.  The  incident  plane  wave  travels  from  left  to 
right  parallel  to  the  horizontal  axis.  Note  that  the  vertical  and  the  horizontal 
scales  have  been  normalized  by  the  wavelength. 

Interpretation  of  plots  labeled  (b)  requires  some  care.  To  explain  what 
different  segments  of  these  plots  stand  for,  we  define  l'nm  as  follows, 

°sfra[Tn)  u,(Tn) 

°/ra|  rj  o.frjTj  u,(7n) 

where  7a  is  a  point  on  the  receiver  line.  Therefore,  Unn  is  the  singly  scattered 
wave  from  object  n,  and  l:nm  is  the  doubly  scattered  wave  generated  by  the 
incident  field  first  striking  the  object  n,  the  resulting  scattered  field  then 
hitting  object  m,  and  the  subsequently  scattered  field  reaching  the  receiver  line. 
Plots  labeled  (b)  show  the  magnitudes  of  lTu  ,  lTjo  ,  IN,  and  l:o2  over  the 
receiver  line  (the  vertical  line  on  the  right  in  plots  a),  simply  plotted  one  after 
another  on  the  same  axis.  In  plots  labeled  (c)  the  solid  line  shows  the 

o  o 

magnitude  of  u3d  =  V  l'nm  (the  singly  plus  doubly  scattered  field) 

n  =1  m  =  l 

O 

whereas  the  dashed  line  represents  the  magnitude  of  u3  =  V  Unn  (the  singly 

n  =  1 

scattered  field).  In  plots  labeled  (d),  the  solid  line  represents  the  phase  of  u.d, 
whereas  the  dashed  line  represents  the  phase  of  us.  The  phase  plot  is  shown 
only  over  the  middle  part  of  the  receiver  line  in  order  to  avoid  the  oscillations 
near  the  ends  of  this  line  which  detract  from  the  effects  we  wish  to  highlight. 
It  should  be  noted  that  the  subscript  I  is  always  used  for  the  cylinder  residing 
on  the  left  or  at  the  bottom  (compared  to  the  other  one).  Also  the  starting 
point  of  horizontal  axes  in  plots  labeled  (b)  and  (c)  corresponds  to  the  top  most 
point  on  the  receiver  line  depicted  in  (a). 

The  cases  that  we  have  studied  can  be  divided  into  two  categories.  In  one 
the  cylinders  are  in  line  with  the  direction  of  illumination,  and  therefore  are 
blocking  each  other.  In  the  other  category,  the  cylinders  are  not  in  the  shadow 
of  their  partners.  All  the  studies  in  diffraction  imaging  assume  either  the  Born 
or  the  Rytov  approximations  as  the  first  st  op.  Although  the  Rytov 
approximation  is  considered  to  be  superior  to  the  Born  approximation  to  some 
extent,  they  both  impose  a  severe  restriction  on  the  maximum  allowable 
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change  of  refractive  index  of  the  object.  The  Born  approximation  implies  th<it 
the  total  field  inside  the  object  can  be  approximated  by  the  incident  field 
However,  the  doubly  scattered  field  U12  (the  scattered  field  of  the  second 
cylinder  when  the  incident  field  is  assumed  to  be  the  scattered  field  of  the  other 
cylinder)  will  be  negligible  only  if  the  total  field  incident  on  the  second  cylinder 
is  close  to  the  original  incident  field.  Therefore  we  expect  that  long  before  the 
Born  approximation  breaks  down  for  each  cylinder,  the  effect  of  double 
scattering  should  become  large  enough  to  distort  the  received  field.  This  is 
borne  out  by  the  computer  simulations  discussed  in  what  follows. 

Figs.  2.15,  2.16  and  2.17  show  the  results  for  the  refractive  index  changes 
of  1,  3  and  5  percent,  respectively.  The  cylinder  diameters  were  6X  for  all 
these  cases.  In  plots  labeled  (b),  the  reader  may  wonder  about  the  large 
magnitude  of  U12  in  relation  to  that  of  lT,,.  Since  the  refractive  index 
variations  in  the  object  are  small,  the  reader  might  rightfully  expect  the 
magnitude  of  the  summation  U,,  + 1  )2  t0  be  close  to  that  of  Ujj,  and  might 
therefore  find  a  contradiction  in  the  large  values  of  U,o  shown.  However  there 
is  no  contradiction.  Since  the  summation  of  Un  and  U12  is  in  the  complex 
domain,  even  if  (he  magnitude  of  U12  is  larger  than  lrn,  the  magnitude  of  the 
summation  of  the  two  may  indeed  be  very  close  to  that  of  Un.  This  is 
illustrated  in  Fig.  2.18,  where  we  have  assumed  a  large  angle,  0,  between  l',, 
and  l ’|2.  This  angle  can  be  expected  to  be  large  when  dealing  with  object 
components  that  are  large  compared  to  a  wavelength,  this  being  caused 
primarily  bv  the  differences  in  the  speed  of  propagation  between  the  second 
component  and  the  background.  To  show  that  indeed  this  is  the  case,  we  have 
chosen  Fig.  2.17,  which  is  the  worst  example  amongst  Figs.  2.15,  2.16  and  2.17 
in  the  sense  that  it  is  characterized  by  the  largest  |  U12|  /|  U,,|  ,  and  plotted  in 
Fig.  2.19  the  magnitude  of  the  fields  U,,  and  l’u  +U,2  on  the  same  scale.  As 
was  expected,  the  difference  between  |  Un+U12|  and  |  Un|  is  small. 

In  Figs.  2.15,  2.16  and  2.17,  a  simple  comparison  between  the  solid  and 
dashed  lines  (single  plus  double  scattered  field  versus  the  single  scattered  field) 
in  plots  labeled  (c)  and  (d)  reveals  that  a  3  percent  change  in  refractive  index  is 
enough  for  the  double  scattering  effect  to  change  the  total  scattered  field 
completely.  Although  the  magnitude  of  the  field  does  not  change  considerably, 

the  deviation  in  phase  is  as  high  as  -y.  Usually  the  Born  approximation  is 

assumed  to  hold  when  the  change  in  the  phase  of  the  incident  field  going 
through  the  object  is  less  than  2rr.  In  figures  2.9,  2.10  and  2  11  the  approximate 
phase  shift  for  each  cylinder  at  a  point,  on  the  axis  joining  the  center  of  the 


n  :  I  01 


•o*: r ia«  or*  «ci!yc*  uni  *mm»*  or*  «ct:v«i  uni 

(c)  (d) 


Figure  2.15 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second- 
order  contributions  are  included  in  one  case  and  neglected  in  the  other. 
Diameters  and  refractive  indices  of  the  cylinders  are  6  X  and  1.01,  respec¬ 
tively.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  fM,  C12,  f.'21 
and  U 22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.16 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  6  X  and  1.03,  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  £/,,,  Ul2,  U2l  and 
0^2  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.17 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  6  X  and  1.05,  respective¬ 
ly.  a)  Geometrv  of  the  simulation,  b)  Magnitudes  of  L’n,  br12,  br2i  and 
components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.18 

A  special  case  of  the  summation  of  different  order  vectors  of  scattered 
fields.  Even  if  the  magnitude  of  the  second-order  contribution  (f/,2)  is 
large  compared  with  that  of  (£/n),  the  magnitude  of  the  summation 
( (/ j i  +  (/ 12)  may  still  be  close  to  that  of  (f/,,),  especially  when  the  phase 
angle  6  is  large. 
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Figure  2.19 

The  magnitude  of  T',,  in  Figure  9.17b  is  represented  by  the  dashed  line, 
while  the  solid  line  shows  the  magnitude  of  the  summation  of  Un  and 
Vn. 
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cylinders  is  respectively  ,  —■  and  Therefore  the  total  scattered  field 

gets  '’distorted”  long  before  the  phase  shift  along  the  length  of  the  objects 
reaches  the  limit  of  2t.  In  Figs.  2.20,  2.21  and  2.22  the  change  in  refractive 
index  is  one,  two  and  three  percent,  but  the  objects  are  taken  to  be  larger  in 
size,  9X  to  be  precise.  This  time  a  two  percent  change  seems  to  be  the  upper 

limit.  The  phase  shift  in  each  cylinder  is  now  -7-,  —  and  -77-. 

3  3  3 

In  Figs.  2.23  and  2.24  we  have  changed  the  position  of  the  cylinders  to  see 
the  effect  when  they  are  not  in  each  other  s  shadow.  The  refractive  index  of 
the  cylinders  is  taken  to  be  1.05  and  2.0,  respectively.  It  is  quite  clear  that 
multiple  scattering  has  a  very  small  effect  on  the  result  in  this  case,  even  for 
very  high  refractive  indices  of  the  object  components.  This  is  expected  because 
the  object  components  are  large  and  thus  forward  scattering.  However,  even 
when  the  size  of  the  cvlinders  is  decreased  in  Figs.  2.23  and  2.26,  no  significant 
change  is  observed  (refractive  indices  were  respectively,  1.05  and  2.0  in  those 
figures).  Since  the  scatterers  are  not  directional  in  this  case,  the  scattered 
energy  spreads  out  in  all  directions  every  time  a  wavefront  hits  a  component. 
The  multiple  scattering  effects  are  minimal  in  this  case  because  the  energy  jn 
the  doubly  scattered  fields  is  much  lower  than  that  contained  in  the  singly 
scattered  fields,  due  to  the  two  isotropic  scattering  interactions  involved  in  the 
former,  compared  with  a  single  such  interaction  in  the  latter. 

2.3.4  The  New  Procedure  In  an  Attenuating  Medium 

Although  the  algorithms  for  calculating  the  scattered  fields  from  non¬ 
attenuating  cylinders  and  media  were  easily  extended  to  lossy  objects  and 
media,  the  same  is  not  true  for  calculating  the  scattered  fields  from  multiple 
cylinders  in  a  lossy  medium.  This  limitation  is  due  to  the  special  characteristic 
of  the  angular  spectrum  expansion  technique  in  a  lossy  medium,  as  discussed  in 
section  2.1. 

In  the  case  of  an  attenuating  medium,  it  is  not  possible  to  associate  the 
Fourier  transform  of  the  field  with  uniTorm  attenuating  plane  waves.  Tim 
Fourier  components  represent  non-uniform  plane  waves  which  possess  different 
directions  of  propagation  and  attenuation.  The  interaction  of  these  non- 
uniform  plane  waves  with  cylindrical  objects  can  not  be  formulated  as  before. 

Assume  the  scattered  field  from  a  cylinder  is  calculated  on  the  line  x  =0 
and  is  expanded  into  non-uniform  plane  waves  propagating  toward  a  second 
cylinder.  The  aim  is  to  formulate  the  interaction  of  these  non-uniform  plane 
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Figure  2.20 

A  comparison  between  the  Selds  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  9  X  and  1.01,  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  l'u,  t',0,  Uy  and 
V 22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.21 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  9  X  and  1.02,  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  frn,  L  j2,  b;o,  and 
U 22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  e,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  wheD  these  contributions  are  neglected. 
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Figure  2.22 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  9  X  and  1.03,  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  £/,,,  t/12,  (/21  and 
U22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.23 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  6  \  and  1.05,  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  t'.l-  ^  t/21  and 
U 22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.24 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  6  X  and  2.0.  respective¬ 
ly.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  (/,,,  U{n,  t/2.  and 
U 22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.25 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  0.4  X  and  1.05,  respec¬ 
tively.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  Uu,  (V12,  (/„, 
and  i  ;2  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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Figure  2.26 

A  comparison  between  the  fields  on  the  receiver  line  when  the  second  ord¬ 
er  contributions  are  included  in  one  case  and  neglected  in  the  other.  Di¬ 
ameters  and  refractive  indices  of  the  cylinders  are  0.4  X  and  2.0,  respec¬ 
tively.  a)  Geometry  of  the  simulation,  b)  Magnitudes  of  L'u,  f/12,  i'2l 
and  L'22  components  of  the  scattered  field  on  the  receiver  line,  plotted  one 
after  another.  c,d)  Magnitude  and  phase,  respectively,  of  the  scattered 
fields  are  shown  by  the  solid  line  when  the  second  order  contributions  are 
included,  and  by  the  dashed  line  when  these  contributions  are  neglected. 
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waves  with  the  second  cylinder, 
plane  waves  can  be  written  as 

,  i  _  r  -  j  (  k.x  +  k,v 

u,(x.y)  -  1'0  eJ 

where  kv  has  a  real  value.  kx 
relationship 


The  mathematical  expression  for  one  of  these 

C2.51) 

would  be  complex  because  of  the  following 


K 


-  Kr  +  i k 


If  the  position  of  the  point  (kx  ,kv)  is  defined  by  (k,.0)  in  cylindrical  coordinates 
Eq.  (2.51)  can  be  rewritten  in  the  following  form 


jk,rco«(0~0)  j  (j  kx,  r  )  oo?0 


u ifr)  =  l'oeJ  ■•'■-"•  e'  '  (2.53) 

where  (r.o)  are  the  coordinates  of  the  point  (x.y)  in  cylindrical  coordinates. 

To  solve  the  boundary  value  problem,  one  has  to  expand  the  incident  fit 
in  terms  of  cylindrical  basis  functions.  The  result  is 


(2.54) 


’sc.  ft 

ffifr)  =  <-;o 

v;  £n  jn  COS  n(0-o)l  Jn(k,r) 

X 

v;  (n  C os(nd)  In(kx  r) 

n  =  0  v  ’ 

n  =0 

where  In  (— jnJn)  is  the  modified  Bessel  function. 


Because  of  the  complex  nature  of  the  expression  in  Eq.  (2.54),  it  is  quite 
difficult,  if  not  'mpossible,  to  impose  the  boundary  conditions  on  Eqs.  (2.16) 
and  (2.54)  and  solve  for  the  coefficients  An  and  Bn.  Therefore,  we  can  not 
account  for  the  interaction  of  the  plane  waves  generated  by  the  angular 
spectrum  expansion  with  the  object  components.  This  conclusion  suggests  that 
the  extension  of  our  algorithm  for  the  calculation  of  scattered  fields  from 
multi-component  objects  to  the  case  of  a  lossy  medium  is  not  feasible. 


2.4  Test  of  Numerical  Results 

A  simple  test  of  validity  or  the  simulation  results  is  to  check  whether  the 
calculated  scattered  fields  satisfy  the  wave  equation.  Due  to  the  Huygen  s 
principle  when  a  field  is  known  over  a  line  in  a  two  dimensional  space,  it  is 
known  over  the  whole  space.  However,  to  test  a  numerical  technique,  one  has 
to  check  the  validity  of  the  result  over  the  range  of  interest  and  prove  the 
numerical  stability  of  the  technique.  To  test  our  computer  programs  for 
calculation  of  scattered  fields,  we  made  use  of  the  formulation  of  propagation 
of  the  angular  spectrum  (see  section  2.1).  We  first  calculate  the  scattered  fields 
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along  the  line  x^Xj  and  then  propagate  it  to  the  line  x  =Xo,  using  the 
formulation  in  Eq.  (2.9).  This  field  is  then  compared  to  the  one  directly 
calculated  by  our  program  along  the  line  x=x2.  We  choose  the  two  ranges  x, 
and  x2  to  be  at  the  two  extremes  of  the  range  of  interest,  to  make  sure  that  the 
field  generated  by  our  program  in  that  region  satisfies  the  wave  equation. 
Although  there  is  a  possibility  that  the  fields  calculated  do  satisfy  the  wave 
equation  but  do  not  correspond  to  the  exact  problem  the  algorithms  are 
designed  for,  one  should  consider  its  very  low  probability  . 

The  solid  line  in  Fig.  2.27(a)  shows  the  magnitude  of  the  scattered  field  as 
calculated  by  our  program  along  a  line  15  wavelengths  from  the  center  of  a 
cylinder  when  it  is  illuminated  by  a  plane  wave.  The  dashed  line  shows  the 
field  which  is  first  calculated  4  wavelengths  from  the  cylinder  and  then 
propagated  to  a  distance  of  15  wavelengths  from  the  object  using  the 
formulation  in  Eq.  (2.9).  The  phase  information  is  reflected  in  Fig.  2.27(b). 
The  attenuation  in  the  medium  is  assumed  to  be  3.5  db/cm.  Notice  that  the 
two  calculated  fields  match  exactly  and  the  dashed  line  may  not  even  be 
recognized  in  this  figure. 

In  Fig.  2.28,  we  have  shown  the  same  information  as  in  Fig.  2.27,  except 
that  the  incident  plane  wave  has  been  replaced  by  a  point  source.  The  point 
source  was  positioned  on  a  line  parallel  to  the  receiver  line  10  wavelengths 
away  from  the  cylinder  and  3  wavelengths  off  axis.  Again  the  two  calculated 
fields  are  very  close. 

For  testing  the  program  that  calculates  the  scattered  field  from  a  multi- 
component  object,  we  chose  the  geometry  in  Figs.  2.15,  where  the  components 
block  each  other,  and  the  one  in  2.23,  where  no  blocking  occurs.  We  chose  the 
incident  field  to  be  a  plane  wave;  however,  a  point  source  illumination  could 
have  been  used  as  well.  In  Figs.  2.29  and  2.30  we  have  shown  the  same 
information  as  plotted  in  Figs.  2.27  and  2.28,  while  using  the  geometry  of  Fig. 
2.15,  with  the  refractive  indices  of  the  components  being  1.1  and  1.5, 
respectively.  Figs.  2.31  and  2.32  are  generated  for  the  geometry  used  in  Fig. 
2.23. 

Since  the  fields  calculated  by  the  method  of  propagation  of  the  angular 
spectrum  are  very  close  to  the  ones  calculated  directly,  we  have  been  persuaded 
that  our  computer  programs  are  error-free. 
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Figure  2.27 

Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  a  cylinder  il¬ 
luminated  by  a  plane  wave,  measured  on  the  receiver  line  at  a  distance  of 
15  X  from  the  cylinder.  Solid  line  represents  the  magnitude  of  the  field 
when  it  is  calculated  directly,  while  the  dashed  line  corresponds  to  the  case 
where  the  field  is  first  calculated  at  a  distance  of  5  X  from  the  cylinder  and 
then  propagated  forward  to  the  receiver  line.  Diameter,  refractive  index 
and  attenuation  constant  of  the  cylinder  are  6  X,  1.1,  and  4.0  db/cm, 
respectively.  The  attenuation  constant  of  the  medium  is  3.5  db/cm. 
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Figure  2.28 

Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  a  cylinder  il¬ 
luminated  by  a  point  source  3  X  off  axis  and  10  X  away  from  the  cylinder, 
measured  on  the  receiver  line  at  a  distance  of  15  X  from  the  cylinder. 
Solid  line  represents  the  magnitude  of  the  field  when  it  is  calculated  direct¬ 
ly,  while  the  dashed  line  corresponds  to  the  eas_  where  the  field  is  first  cal¬ 
culated  at  a  distance  of  5  X  from  the  cylinder  and  then  propagated  for¬ 
ward  to  the  receiver  line.  Diameter,  refractive  index  and  attenuation  con¬ 
stant  of  the  cylinder  are  6  X,  1.1,  and  4.0  db/'em,  respectively.  The  at¬ 
tenuation  constant  of  the  medium  is  3.5  db/cm. 
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Figure  2  20 

Magnitude  (a)  and  phase  { b )  of  the  scattered  field  from  the  two-cylinder 
object  shown  in  Figure  9.15a  measured  on  the  receiver  line  at  a  distance  of 
20  X  from  the  origin.  Solid  line  represents  the  magnitude  of  the  field  when 
it  is  calculated  directly,  while  the  dashed  line  corresponds  to  the  case 
where  the  field  is  first  calculated  at  a  distance  of  10  X  from  the  origin  and 
then  propagated  forward  to  the  receiver  line  The  refractive  index  of  the 
cvlinders  is  1  1 
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Figure  2.30 

Magnitude  (a)  and  phase  (b)  of  the  scattered  Geld  from  the  two-cylinder 
object  shown  in  Figure  9  15a  measured  on  the  receiver  line  at  a  distance  of 
20  X  from  the  origin  Solid  line  represents  the  magnitude  of  the  Geld  when 
it  is  calculated  directly,  while  the  dashed  line  corresponds  to  the  case 
where  the  Geld  is  first  calculated  at  a  distance  of  10  X  from  the  origin  and 
then  propagated  forward  to  the  receiver  line.  The  refractive  index  of  the 
cylinders  is  1.5. 
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Figure  2.32 

Magnitude  (a)  and  phase  (b)  of  the  scattered  field  from  the  two-cylinder 
object  shown  in  Figure  9.21a  measured  on  the  receiver  line  at  a  distance  of 
20  X  from  the  origin.  Solid  line  represents  the  magnitude  of  the  field  when 
it  is  calculated  directly,  while  the  dashed  line  corresponds  to  the  case 
where  the  field  is  first  calculated  at  a  distance  of  10  X  from  the  origin  and 
then  propagated  forward  to  tbe  receiver  line.  The  refractive  index  of  the 
cylinders  is  1.5. 
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CHAPTER  3 
INVERSE  PROBLEM 


In  general  inverse  scattering  is  a  much  more  complicated  problem  than  the 
forward  problem.  Inverse  scattering  techniques  are  usually  accompanied  by 
severe  restrictions  on  the  different  parameters  of  the  system  under 
consideration. 

Small  perturbation  solutions  for  diffraction  reconstruction  techniques  were 
initiated  in  1969  by  Wolf  [Wolf69].  He  proved  that  when  a  monochromatic 
plane  wave  is  incident  on  an  object,  the  Fourier  transform  of  the  scattered  data 
measured  on  a  receiving  planar  aperture  and  the  Fourier  transform  oF  the 
object  function  are  linearly  related.  Dandliker  and  Weiss  obtained  similar 
results  for  three-dimensional  reconstruction  [Dandliker70],  Carter  used  this 
method  for  reconstructing  refractive  index  distribution  by  using  holograms  to 
record  the  scattered  fields  in  the  Fresnel  zone  of  a  semi-transparent  object 
[Carter70],  In  spite  of  many  experimental  problems,  his  results  were  quite 
promising.  Iwata  and  Nagata  addressed  some  of  the  limitations  of  this 
approach  and  proposed  a  cylindrical  receiver  line  vis-a-vis  the  straight  line  (as 
suggested  by  Wolf)  for  small  objects  [Iwata75], 

Some  theoretical  inverse  procedures  for  reconstructing  weakly 
inhomogeneous  objects  have  been  recently  derived  by  several  researchers.  A 
three-dimensional  inverse  scattering  scheme  for  a  spherical  recording  aperture 
enclosing  the  object  was  derived  by  Ball  et  al.  [BallSO].  They  considered 
monochromatic  spherical  wave  sources  for  the  object  insonification  and  point 
receivers  for  the  measurement  of  the  scattered  field.  Norton  and  Linzer 
[Norton81]  have  considered  broadband  insonification  to  derive  reconstruction 
methods  for  planar,  spherical  and  cylindrical  geometries. 

Mueller  et  al.  employed  the  method  suggested  by  Wolf  to  implement  a 
cross-sectional  reconstruction  scheme  for  imaging  the  refractive  index  of  objects 
[Mueller79j.  They  have  shown  that,  a  coverage  of  the  object's  Fourier 
transform  function  can  be  achieved  by  measuring  the  scattered  field  on  a 
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receiving  line  for  different  rotational  positions  of  the  object  over  360  degrees. 
This  method  provides  the  values  of  the  object’s  Fourier  transform  on  a  set  of 
arcs,  while  for  reconstructing  the  object  function  through  the  discrete  Fourier 
transformation  on  a  digital  computer  the  values  are  required  on  a  two- 
dimensional  rectangular  grid.  Different  interpolation  techniques  have  been  used 
for  this  purpose  [Pan83,  Kaveh79a].  Devaney  [Daveaney82]  has  derived  a 
method,  which  he  has  called  the  Fourier  back-propagation  algorithm  (similar  in 
spirit  to  Fourier  Back-Projection  for  straight  ray  tomography),  to  avoid 
interpolation  errors.  A  drawback  of  this  approach  is  its  extensive  demand  on 
computer  time.  A  similar  approach  has  been  followed  independently  by 
Soumekh  et  al.  [Soumekh83j. 

Nahamoo  and  Kak  (Nahamoo82,  Nahamoo84j  generalized  the  technique  of 
diffraction  tomography  as  presented  in  [Mueller79]  to  derive  a  new  scheme 
which,  in  principle,  permits  the  use  of  any  type  of  insonification,  as  opposed  to 
only  the  plane  wave  type  which  is  experimentally  difficult  to  generate.  With 
this  approach,  the  requirement  that  the  object  be  viewed  over  360  degrees  is 
replaced  by  the  need  for  merely  two  views,  which  ideally  should  be  90  degrees 
apart  (however,  that  is  not  strictly  necessary).  For  each  of  these  two  views  the 
method  requires  independent  traversal  movements  for  both  the  transmitter  and 
the  receiver  on  two  parallel  lines  in  the  plane  of  the  cross  section. 

Devaney  has  recently  taken  an  approach  similar  to  that  of  Nahamoo  and 
Kak  in  applications  dealing  with  geophysical  imaging  [Devaney84].  He  has 
considered  the  case  where  the  transmitter  and  receiver  arrays  are  perpendicular 
to  each  other  and  has  derived  the  corresponding  point  spread  functions. 

Kenue  and  Greenleaf  [Kenue82j  extended  the  approach  of  Mueller  et  al. 
[Mueller79]  to  the  multi-frequency  case.  They  have  shown  that  by  measuring 
the  projections  at  more  than  one  frequency,  one  can  reduce  the  number  of 
views  required  and  thereby  higher  speeds  in  data  collection  may  be  achieved. 
Note  that  the  total  amount  of  data  collected  remains  basically  the  same, 
because  for  each  rotational  position  of  the  object,  each  projection  must  now  be 
recorded  at  different  frequencies.  It  is  suspected  that  in  the  experimental 
implementation,  the  frequency  dependence  of  tissue  attenuation  might  become 
a  source  of  difficulty  for  this  method. 

Recently  Johnson  et  al.  [Johnson83,  Tracy83,  Johnson81)  have  suggested  a 
method  based  on  solving  systems  of  nonlinear  algebraic  equations  that  are 
derived  by  applying  the  method  of  moments  to  a  bine  basis  function  expansion 
of  the  fields  and  scattering  potential.  Although  it  seems  that  the  limitations  of 
this  approach  is  quite  similar  to  that  of  the  Born  approximation,  the  domain  of 
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applicability  of  this  approach  still  remains  to  be  proved. 

In  section  3.1,  we  have  reviewed  the  bascis  of  the  interaction  of  a  wave 
with  tissue,  the  integral  representation  of  the  scattered  field  and  the  two  well 
known  weak  assumptions:  the  Born  and  the  Rytov.  The  conventional 
Diffraction  Reconstruction  technique  is  discussed  in  section  3.2,  and  the 
extension  of  the  technique  to  the  case  of  an  attenuating  medium  examined. 
We  review'  the  mathematical  derivation  of  the  Synthetic  Aperture  technique  in 
transmission  mode  [Nahamoo82]  in  section  3.3.  The  theory  is  extended  to  the 
reflection  mode  in  that  section  and  it  is  proved  that  if  the  reconstruction 
technique  is  utilized  in  both  the  transmission  and  the  reflection  modes,  for  a 
single  view  of  the  object  one  can  recover  the  object’s  Fourier  transform  inside  a 
circle  of  radius  2k0  (k0  is  the  wavenumber).  Finally  in  section  3.3.3  the  effect 
of  attenuation  on  this  technique  is  considered. 

3.1  Formulation  of  the  Problem 

The  interaction  of  sound  with  a  medium  forces  the  particles  of  that 
medium  to  be  displaced  from  their  equilibrium  positions.  The  elastic  restoring 
forces  combined  with  the  inertia  of  the  particles  lead  to  an  oscillatory  motion 
of  the  medium.  A  mathematical  description  for  the  propagation  of  sound  in  a 
medium  can  be  derived  by  formulating  the  principles  of  the  elastic  properties  of 
the  medium  and  setting  up  an  equation  of  motion  for  the  particles  in  the 
medium. 


3.1.1  Interaction  of  Wave  with  Tissue 


We  only  consider  the  longitudinal  component  of  the  sound  wave,  assuming 
the  traversal  components  get  damped  because  of  the  fluid  structure  of  the 
medium  and  the  object  under  consideration  This  is  usually  true  for  tissue 
because  of  its  high  fractional  water  content.  Based  on  this  model,  the  equation 
of  motion  is  given  by  [Morse68], 

pm  =  -vpfr.t)  .  (3.i) 

at 

where  V  and  P  are  velocity  and  pressure  fields,  respectively  and  p  is  the  local 
density  of  the  fluid.  The  equation  of  state  of  the  fluid  near  the  equilibrium 
state  (small  amplitude  disturbances)  is 


-  pff.t)  Kfr.t) 


where  k  is  the  complex  compressibility  of  the  viscous  fluid.  The  local  variation 
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of  density  is  related  to  the  flow  of  fluid  at  each  point  by  the  equation  of 
continuity  as  given  by 

MlA  =  -  pfrt)  V  ■  V(7, t)  (3  3) 

dt 

By  combining  the  three  equations  of  motion  (Eq.  (3.1)),  state  (Eq.  (3.2)),  and 
continuity  (Eq.  (3.3)),  we  come  up  with  the  following  wave  equation. 


pm 


VPff.t)  J  =  /c(f,t) 

dpm. 

-(-  dt 

pm 


pm  *(F,t) 


pm 


The  hist  term  on  the  right  hand  side  of  this  equation  is  the  difference  between 
the  percentage  change  in  time  of  p  k  and  p  multiplied  by  the  percentage 
change  in  time  of  p.  This  term  being  small  compared  to  the  first  term  on  the 
right  hand  side  of  Eq.  (3.4)  is  usually  neglected.  For  a  single  frequency  case 
with  P(T,t)  =  u(T)  e^1,  we  obtain 

V  •  —  Vu(7)  +  k(7)  u)'“  ufr)  =  0  (3.5) 

p(n 

The  wavenumber  is  defined  as: 

kco  “  Po  Ko  'J  »  (3-6) 

where  p0  and  k0  are  the  density  and  the  compressibility  of  the  fluid  at 
equilibrium,  respectively.  When  the  fluid  is  lossy  the  compressibility,  k0, 
becomes  complex  and  so  does  kco. 

If  we  introduce  the  notations 


- 


ptf)  ~  Po 


7,-m  = 


K’cn  -  k0 


(3.7b) 


and  make  use  of  definition  of  kro  in  Eq.  (3.6),  we  can  rewrite  Eq.  (3.5)  in  the 
following  form 


(3.8) 


V  '  (  (  1  -  7,(T)  ]  Vuprl  )  +  k„2  |  1  +  7„(J)  ]  <i(T)  =  0  . 

By  massaging  the  equation  above,  we  get  the  final  form  of  the  wave  equation 

(V2  +  IO  u(F)  =  -  k„2[n2(r)-l|  U(r)  +  -e0-  VV(T)  ■  Vufr)  ,  (3. 

Po 


where  the  refractive  index  n(7)  is  defined  as 


n(7)  - 


(3.10) 


In  most  practical  situations  the  second  term  on  the  right  hand  side  of  (3.9) 
is  neglected  and  the  following  simpler  linear  inhomogeneous  Helmholtz  equation 
is  used 


(V2  +  kcc2)u(T)  =  -F(7)  ufr)  , 
where  the  source  function  Ffr)  is  defined  as 
F(f)  =  kco2  [n2(T)  -  1] 


(3.11) 


(3.12) 


In  all  the  simulations  considered  in  our  study,  the  density  function  p(r)  is 
assumed  to  be  equal  to  p0,  so  that  7  (7)  becomes  zero  and  the  last  term  on  the 
right  hand  side  of  Eq.  (3.9)  vanishes. 


3.1.2  An  Integral  Solution 

The  aim  is  to  convert  the  differential  equation  (3.11)  to  an  integral 
equation.  An  elegant  description  of  this  problem  can  be  found  in  [Morse53], 
while  in  here  we  only  bring  a  brief  summary  of  the  approach  similar  to  that 
covered  in  (Nahamoo82|. 

The  wave  equation  as  given  in  Eq.  (3.11)  does  not  provide  a  complete 
definition  of  the  problem  unless  it  is  accompanied  by  the  boundary  conditions 
which  specify  the  behavior  of  the  field  on  the  surface  that  encloses  the  medium. 
The  boundary  conditions  can  be  expressed  in  terms  of  the  pressure  field  and 
the  normal  velocity  at  the  boundary.  Since  the  velocity  is  related  to  the 
derivative  of  the  pressure  according  to  Eq.  (3.1),  the  pressure  field  and  its 
normal  derivative  at  the  boundary  are  usually  employed  to  express  the 
boundary  conditions. 

In  practical  cases,  the  boundary  consists  of  partially  rigid  or  flexible 
materials,  while  on  some  part  of  it  a  velocity  distribution  is  generated  by  the 
transmitting  transducer.  A  general  representation  of  such  boundary  conditions 
is  given  by 


^B^b)  '  °B^b)  +  zC*b)  uHb)  ~  factive(r*B)  >  (3.13) 

where  7B  is  a  point  on  the  boundary,  7B(7B)  is  the  unit  vector  normal  to  the 
boundary,  z(7B)  represents  the  reaction  impedance  between  the  medium  and  the 
boundary  and  factiveOB)  is  a  function  which  is  non-zero  for  the  active  part  of 
the  boundary.  For  a  transmitting  transducer,  zCie)/jko/)0c0  represents  the 
impedance  of  the  transducer  and  factiveOB)/j^oPoco  *s  ca^e(f  the  driving  velocity 
of  the  transducer. 

The  wave  equation  given  in  Eq.  (3.11)  is  called  inhomogeneous  because  the 
right  hand  side  of  the  equation  is  non-zero  and  represents  a  distributed  source 
of  energy.  The  wave  equation  is  called  homogeneous  if  the  function  F(F)  is 
identically  zero  everywhere.  Since  the  function  factivetrB)  contributes  to  the 
production  of  the  field  almost  in  the  same  manner  as  does  F(r),  it  is 
appropriate  to  name  the  boundary  conditions  similar  to  the  wave  equation. 
Homogeneous  boundary  conditions  are  associated  with  factive(*B)  being 
indentically  zero,  while  a  non-zero  factive(rB)  represents  inhomogeneous 
boundary  conditions. 

The  unique  solution  of  an  inhomogeneous  wave  equation  with  the 
associated  boundary  conditions  is  the  sum  of  two  components.  One  is  the 
solution  of  the  inhomogeneous  equation  for  the  homogeneous  boundary 
conditions,  while  the  other  represents  the  solution  of  the  homogeneous  equation 
for  the  inhomogeneous  boundary  conditions.  Clearly  the  sum  of  these  two 
solutions  will  satisfy  both  the  inhomogeneous  wave  equation  and  the 
inhomogeneous  boundary  conditions. 

The  solution  for  the  homogeneous  wave  equation  under  inhomogeneous 
boundary  conditions  is  the  field  generated  by  the  transmitting  transducer, 
which  we  call  Uj.  The  solution  for  the  inhomogeneous  wave  equation  under 
homogeneous  boundary  conditions  is  obtained  by  adding  the  contribution  of  all 
the  differential  elements  of  the  source  function,  F.  The  field  at  point  7  caused 
by  a  point  source  at  7Q  is  the  well-known  free  space  Green’s  function  g(7[70) 
which  satisfies  the  following  wave  equation 

(V2  +  kj)  g(7j  7,)  =  -<5(7-70)  (3.14) 

The  field  generated  by  the  distributed  source  function  F  is  the  integral  of 
F(70)  gfrj70)  over  the  scattering  object  Therefore,  the  integral  representation 
of  the  differential  equation  in  Eq.  (3.11).  can  be  written  as 
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us(r)  =  /  F(T0)  u(f0)  g(r]T0)  dr0  , 

(3.15) 

object 

where  u3  is  the  scattered  field  defined  by 

us(r)  =  u(F)  -  ujfr) 

(3.16) 

Since  the  total  field,  u,  does  appear  inside  the  integral  in  Eq.  (3.15),  the 
solution  is  of  transcendental  form,  which  means  an  inverse  reconstruction 
technique  based  on  this  solution  is  not  feasible.  Some  small  perturbation 
approximation  leads  to  an  expression  for  the  scattered  field  us  as  an  explicit 
function  of  the  inhomogeneities.  Two  of  the  well  known  approximations  are 
the  Born  and  the  Rytov  approximations. 

3.1.3  Born  Approximation 

In  the  first-order  Bom  approximation  the  total  field,  u,  on  the  right  hand 
side  of  Eq.  (3.15)  is  approximated  by  the  incident  field,  u;.  This  is  just;  dable 
only  when  the  object  inhomogeneities  are  small  so  that  the  magnitude  of  the 
scattered  field  inside  the  object  is  negligible  compared  to  that  of  the  incident 
field.  Then  the  solution  is  given  by  the  following  integral 

usffl  =  /  F(70)  Ui(f0)  g(r|70)  dT0  (3.17) 

object 

Second  order  Born  solutions  can  be  obtained  by  replacing  the  total  field, 
u(7)  in  Eq.  (3.15)  by  the  sum  of  the  incident  field  and  the  first-order  scattered 
fields  calculated  using  Eq.  (3.17).  A  repetitive  application  of  this  procedure 
results  in  higher  order  solutions.  It  is  known  that  the  higher  order  solutions 
converge  to  the  exact  solution  if  the  object  inhomogeneities  are  small.  Although 
a  weak  sufficient  condition  for  the  convergence  of  the  solution  exists,  the 
necessary  condition  is  still  unknown  [Hochstadt73j. 

3.1.4  Rytov  Approximation 

An  alternative  to  the  Born  is  the  Rytov  approximation.  In  this  case  the 
incident  and  the  scattered  fields  are  represented  by 

u.fr)  =  e*jn  (3.18a) 

and 

U,(r)  =  e*,!rl  (3.18b) 


Using  these  representations  of  the  fields,  the  corresponding  integral  equation 
can  be  written  as: 

^00  =  -7-  /  (l  V<*s(n)  I  2  +  Ffro)l  u,(f0)  gff]70)  df0  .  (1.19) 

Under  the  first-order  Rytov  approximation,  the  square  of  the  magnitude  of  the 
scattered  field  gradient  is  assumed  to  be  negligible  compared  to  the  object 
function,  so  that  the  term  in  brackets  can  be  approximated  by  Ffr0). 
Therefore  the  final  equation  for  the  scattered  phase  subject  to  the  first-order 
Rytov  approximation  is 

d>sm  =  I  F(70)  U|(70)  g(7|70)  d?0  .  (3.20) 

^iCF)  object 

A  second  order  solution  is  obtained  through  substitution  of  the  gradient  of  the 
first-order  solution  in  Eq.  (3.19).  As  in  the  case  of  higher  order  Born  solutions, 
by  applying  the  same  technique  repetitively  one  can  come  up  with  higher  order 
Rytov  solutions. 

It  is  important  to  note  that,  in  spite  of  the  similarities  between  the  Born 
(Eq.  (3.17))  and  the  Rytov  (Eq.  (3.20))  solutions,  the  approximations  are  quite 
different  [Keller69,  SancerTO].  The  Born  approximation  produces  a  better 
estimate  of  the  scattered  field  amplitude  for  large  deviations  in  the  refractive 
index  of  objects  which  are  small  in  size.  On  the  other  hand  the  Rytov 
approximation  gives  a  more  accurate  estimate  of  the  scattered  phase  for  large 
sized  objects  with  small  deviations  in  refractive  index.  When  the  object  is  small 
and  the  refractive  index  deviates  only  slightly  from  that  of  the  surrounding 
media,  it  is  possible  to  show  that  the  Born  and  Rytov  approximations  produce 
the  same  results  [SlanevSl], 

3.2  Conventional  Diffraction  Tomography 

The  basic  concept  behind  most  of  the  present  Diffraction  reconstruction 
techniques  has  been  derived  by  Wolf  in  1969  [Wolf69],  This  concept,  called  the 
Fourier  Diffraction  Projection  Theorem  by  Pan  and  Kak  [Pan83],  can  be  stated 

as: 

When  an  object  is  illuminated  with  a  plane  wave  as  shown  in  Fig. 

3  1,  the  Fourier  transform  of  th«'  forward  scattered  fields 
measured  on  a  line  perpendicular  to  the  direction  of  propagation 
of  the  wave  (line  TT'  in  Fig.  3.1)  gives  the  values  of  the  two- 
dimensional  Fourier  transform  of  the  object  along  a  circular  arc. 

It  should  be  noted  that  the  Fourier  Diffraction  Projection  Theorem  is  valid 
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only  when  the  object  inhomogeneities  are  weakly  scattering,  or  in  fact  when 
either  the  Born  or  the  Rytov  approximations  hold. 

In  the  following,  we  review  the  derivation  of  this  concept  and  extend  it  to 
the  case  of  an  attenuating  medium,  assuming  the  Born  approximation. 
Consider  the  geometry  shown  in  Fig.  3.2.  We  assume  an  incident  plane  wave 
propagating  along  the  axis  7,  expressed  by 

u,(r)  =  U0eik“Rr,r  ,  (3.21) 

where  f/  is  the  unit  vector  along  the  7  axis  and  R  is  the  rotational 
transformation  applied  to  the  £-7  coordinates  to  position  it  along  the  x-y 
coordinates. 

The  solution  for  the  Green’s  function  in  two-dimensional  space  is  a  Hankel 
function  of  zeroth  order  [Morse53j 
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If  kco  is  complex,  the  Hankel  function  dies  out  exponentially  as  a  function  of 
|7-70|.  This  becomes  obvious  when  one  considers  the  asymptotic  form  of 
this  function. 


Ho0)(z)  -  -ie'2'2’  (3.23) 

z 

This  behavior  agrees  with  our  physical  intuition  of  the  problem,  since  the 
energy  from  a  point  source  in  a  lossy  medium  should  decay  rapidly  as  a 
function  of  the  distance. 

The  Green’s  function  is  now  expanded  in  terms  of  plane  waves  as 
[Morse53], 

:  00  JRXtf-n) 

g(7|?0)  =  -j-  f  -  da  (T~70)  Rri  >  0  ,  (3.24) 

47r  -oo  P 


where 


A  =  (M) 


(3.25) 


P  =  s/Ko2-*2  .  (3.26) 

and  n  and  J  axes  are  in  the  same  directions  as  f  and  r;  are.  The  condition 
(7-70)  >  0  in  Fq.  (3.21)  restricts  the  domain  of  the  plane  wave  expansion 

expressed  by  that  equation  to  7  >  0.  Note  that  when  kro  is  complex,  3  would 
also  be  a  complex  quantity.  1'hus  each  plane  wave  is  attenuated  in  the  7 


direction  as  it  propagates  through  the  half  space  tj  >  0. 

Substituting  Eqs.  (3.21)  and  (3.24)  in  Eq.  (3.17)  and  interchanging  the 
order  of  integrations,  we  have 

usfr)  =  “  /  f  Ff7o)  e^R(X_k^,r<>  d?0da  .  (3.27) 

47r  -00  P  object 

Taking  the  Fourier  transform  of  the  field  on  the  line  r)-r]0}  we  get 

jU  _  \ 

U»  =  ^  F(r  (\-K0f,)\  ,  (3.28) 

where 

CO 

us(a)  =  /  us(^0)  d?  ,  (3.29) 

-00 

and 

+  I  Fff0)e-i?“(?  +  ^dr0  .  (3.30) 

object 

the  ^  and  q[  vectors  are  in  the  direction  of  x  and  y  axes,  respectively. 

For  a  lossless  medium  (real  kco),  F  as  defined  in  Eq.  (3  30)  represents  the 
two-dimensional  Fourier  transform  of  the  object  function  F  (because  ^  would 
be  zero).  Considering  the  relationship  between  a  and  j3  as  given  by  Eq.  (3.26), 
it  is  easy  to  show  that  one  can  recover  the  Fourier  transform  of  the  object 
function  on  the  semicircle  shown  in  Fig.  3.3  by  employing  Eq.  (3.28).  As  the 
object  is  rotated  over  360  degrees,  the  Fourier  domain  is  filled  up  (see  Fig.  3.4), 
and  then  a  reconstruction  can  be  obtained  by  the  Fourier  transform  inversion; 
this  being  the  approach  initially  suggested  by  Mueller  et  al.  for  the 
conventional  Diffraction  Tomography  (Mueller79). 

When  the  medium  is  attenuating,  (J  is  non-zero  in  Eq.  (3.30)  and  that 
equation  defines  a  Laplace  transform  relationship  between  the  functions  F  and 
F.  The  rotation  of  the  object  would  provide  one  with  a  coverage  in  the  four- 
dimensional  Laplace  domain  over  a  two-dimensional  surface.  An  inverse 
Laplace  transformation  is  possible  only  if  cj  is  a  fixed  vector  [Bellman66, 
Krylov69,  Krylov77],  A  close  look  at  Eqs.  (3.25),  (3.26)  and  (3.27)  reveals  the 
fact  that  p  and  cf  are  dependent  variables  and  (j  would  change  whenever  p  is 
altered  Therefore,  the  conventional  reconstruction  technique  can  not  be 
extended  to  the  case  oi  an  attenuating  medium. 
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Figure  3.3 

The  Fourier  transform  of  the  scattered  field  on  the  receiver  line  gives  the 
values  of  the  object’s  Fourier  transform  over  a  semi-circle  in  the  Fourier 
domain. 
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Figure  3.4 

By  rotation  of  the  object  over  360  degrees  the  object’s  Fourier  transform  is 
recovered  over  a  set  of  semi-circles  inside  the  circle  v2 k,. 


The  object’s  cross  section  can  be  recovered  by  taking  the  inverse  Fourier 
transform  of  the  recovered  values  in  the  object  Fourier  domain.  To  use  an 
FFT  based  algorithm  for  the  inverse  transform,  one  has  to  know  the  values  of 
the  function  over  a  square  grid.  This  calls  for  an  interpolation  process  from  the 
cicular  arcs  to  the  square  grid.  The  resulting  interpolation  errors  lead  to 
artifacts  in  the  reconstructed  object.  Different  interpolation  skims  are 
compared  in  [Pan83j.  If  the  number  of  projections  and  the  number  of  sampled 
points  per  projection  are  large  enough,  bilinear  interpolation  seems  to  produce 
relatively  small  artifacts  [Pan83j.  This  can  also  be  justified  by  considering  Fig. 
3.5  for  the  case  where  the  cicular  arcs  and  the  sampled  points  are  very  closely 
spaced. 

It  has  been  shown  by  Devaney  [Devaney83]  that  the  Fourier  Diffraction 
Projection  Theorem  can  also  be  used  to  derive  an  algorithm  that  is 
conceptually  similar  to  the  filtered-backprojection  algorithm  for  x-ray 
tomography.  This  algorithm  eliminates  the  need  for  interpolation  in  the 
frequency  domain  and  he  called  it  the  filtered-backpropagation  algorithm.  In 
the  backpropagation  technique  the  filtered  projections  are  mapped  onto  the 
object  space  via  an  integral  transform  that  is  shown  to  be  the  inverse  of  the 
transform  that  governs  phase  propagation  within  the  Rytov  approximation.  In 
spite  of  the  conceptual  similarities  between  the  backprojection  and  the 
backpropagation  algorithms,  for  its  full  implementation,  for  an  NxN 
reconstruction  with  N  projections,  the  filtered-backprojection  algorithm  needs 
only  N  Fourier  transforms.  To  contrast,  the  filtered-backpropagation 
algorithm  requires  about  N2  Fourier  transforms.  If  the  extensive  demands  on 
computer  time  are  not  a  factor,  one  can,  of  course,  use  the  filtered- 
backpropagation  algorithm.  In  most  practical  applications,  however,  one  may 
have  to  resort  to  procedures  that  are  computationally  more  efficient. 

3.3  Synthetic  Aperture  Diffraction  Tomography 

Nahamoo  and  Kak  [Nahamoo82]  have  derived  a  new  synthetic  aperture 
imaging  technique  that  permits  the  use  of  any  type  of  insonification,  as 
opposed  to  only  the  plane  waves  and  also  requires  only  two  rotational  positions 
of  the  object.  For  each  rotational  position,  the  projection  data  is  collected  by 
independent  traversal  movements  for  the  transmitting  and  the  receiving 
transducers  (see  Fig.  3.6).  Although  ideally  the  angle  between  the  two 
rotational  positions  should  be  90°,  theory  predicts  that  valid  results  should  be 
obtainable,  albeit  with  reduced  spatial  resolution,  even  when  this  condition  is 
not  satisfied.  For  each  rotational  position  of  the  object,  the  data  can  be 
collected  most  efficiently  by  using  arrays  on  both  the  transmit  and  the  receive 
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sides.  The  elements  of  the  transmit  array  could  be  fired  sequentially  and  for 
each  such  firing  the  received  field  would  be  measured  over  all  the  elements  of 
the  receiver  array.  In  section  3.3.1,  we  review  the  theory  behind  the  Synthetic 
Aperture  technique  in  transmission  mode,  as  proposed  by  Nahamoo  and  Kak. 
Later  in  section  3.3.2,  we  prove  that  if  the  system  is  utilized  in  both  the 
transmission  and  the  reflection  modes,  only  one  rotational  position  of  the  object 
will  be  sufficient  for  the  reconstruction  of  the  object's  cross  section.  We 
examine  the  Synthetic  Aperture  technique  in  the  presence  of  attenuation  in  the 
last  section  (3.3.3). 

3.3.1  Transmission  Mode 

Consider  the  geometry  in  Fig.  3.6.  When  the  transmitter  is  located  at  the 
point  (xt,yt),  let  the  incident  field  U;  and  the  scattered  field  u3  at  a  point 
T  =  (x,y)  be  represented  by  Uj(7;yt)  and  us(7;yt)  respectively.  The  incident  field 
Uj(T;0)  can  be  represented  by  the  angular  spectrum  expansion  (see  section  2.1). 


ui(T;°)  =  /  At(ky;0)  ejR  (7"x‘i)  dk, 


(3.31) 


where  x  is  the  unit  vector  along  the  x  axis.  The  wave  vector  R  is  given  by  (see 
Fig.  3.7  for  real  kco) 

fv  =  (kx,ky)  (3.32) 


kx  =  V 1 kco*  "  V 


(3.33) 


The  function  At(kx;0)  represents  the  complex  amplitude  of  the  angular 
spectrum  of  the  incident  field  at  the  line  x=xt  when  the  transmitter  is 
positioned  at  y=0.  This  function  is  the  Fourier  transform  of  the  incident  field 
on  the  transmit  line  as  given  by 


At(k  ;0)  -  /  Uj(xt,y;0)  e"jk,y  dy 


(3.34) 


The  angular  spectrum  of  the  field  when  the  transmitter  is  located  at  y=yt  is 
d ('fined  as 


At(kryt)  -  /  u i ( x t . y : y t )  e~jk'y  dy 


Since  ui(xl.y-yt;yl)  is  the  same  as  U;(xt.y,0),  we  have 


(3.35) 
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At( ky ;y t)  =  At(ky;0)  e  jk^  .  (3  36) 

Finally  the  following  expansion  for  the  field  is  found.  Now  the  field  at  point  7 
when  the  transmitter  is  at  (xt,yt)  can  be  expressed  as 


UitfyJ  =  “  /  At(ky)  e“Jlc'y‘  e- 


jkyy» 


where  we  have  replaced  At(ky;0)  by  At(ky)  for  convenience  in  presentation. 

Going  back  to  the  basic  integral  equation  in  Eq.  (3.17),  the  Green’s 
function  is  represented  by  an  angular  spectrum  expansion.  The  expansion  for 
x  >xQ  is  given  by 


;  *  JA-(7-r.J 

^l7"1  =  £  /  -T" ' 

-  X) 

where  the  wave  vector  A  is  given  by  (see  Fig.  3.7  for  real  kco) 
A  =  (a,/?) 


(3.38) 


(3.39) 


a  -  \Ae?  “  ^ 


(3.40) 


Let  us(y r:y t)  represent  the  scattered  field  at  the  point  (xr,yr)  on  the  receiver  line 
when  the  transmitter  is  at  (xt,yt).  We  substitute  equations  (3.37)  and  (3.38)  in 
equation  (3.17)  to  obtain 

i  00  --It  1  JforXr  +  A'r) 

us(yr;yt)  =  f  f(g)  II  At(M e  J  ,y‘  —  dky d^ 

~(“tr)  object  -oo 


fj(A-K)r0  (Jf 


(3-41) 


Switching  the  order  of  integrations  and  rearranging  terms,  we  get 
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AJk  )  F(A-R)  elf9yr  e~jk'y'  d/3  dk 


(3.42) 


where  F  is  the  function  defined  in  (3.30).  Taking  the  Fourier  transform  of  both 
sides  of  this  equation  with  respect  to  yr  and  v,,  we  have 


>T(Tky)  - 
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At( k  )  F(A-K) 


where  ud  d,ky)  has  been  defined  as 


(3.44) 


us(^,ky)  =  If  us(yr;-yt)  e'jk^>  dyr  dyt 

-00 

Let  Ps(xr;xt)  be  the  measured  scattered  field  when  the  transmitter  and  the 
receiver  are  located  at  points  (xt,yt)  and  (xr,yr),  respectively.  We  assume  that 
the  receiver  is  a  flat  transducer  with  a  uniform  sensitivity  over  its  aperture. 
The  effect  of  the  receiver  can  then  be  represented  by  an  aperture  function 
Ir(yr).  The  measured  scattered  field  is  then  given  by 


Ps(yrTt)  -  /  us(y;yt)  Ir(vr~y)  dy 

“OO 

A  Fourier  transform  representation  of  equation  (3.45)  is  given  by 
Ps(/?,ky)  -  A t{8)  us(/?,ky)  , 

where  \(3)  is  the  Fourier  transform  of  the  aperture  function  Ir(yr). 
If  w'e  substitute  equation  (3.43)  in  equation  (3.46),  we  obtain 


(3.45) 


(3.46) 


F(A-K)  = 


-2 jo  ejk»*‘  e'iaXr 
At(ky) 


P.(Akv) 


(3.47) 


Nahammo  and  Kak  [Nahamoo82]  have  shown  that  for  a  lossless  medium, 
for  which  F  is  the  object’s  two-dimensional  Fourier  transform  (see  Eq.  (3.30)), 
Eq.  (3.47)  defines  a  one-to-one  relationship  between  the  data  u3  and  the  object’s 
Fourier  transform,  F,  inside  the  region  shown  in  Fig.  3.8.  The  values  of  the 
object's  Fourier  transform  are  recovered  over  the  semicircle  drawn  by  the  solid 
line  in  Fig.  3.8  centered  at  -lx.  The  coverage  area  shown  in  the  figure  is 
achieved  by  moving  the  center  of  the  solid  semicircle  over  the  dashed  semicircle 
centered  at  the  origin.  Therefore,  for  two  rotational  positions  of  the  object,  90 
degrees  apart,  one  can  recover  the  object’s  Fourier  transform  inside  a  circle  of 
radius  \/2k0,  as  shown  in  Fig.  3.9. 

Similar  to  the  conventional  approach,  the  values  of  the  object’s  Fourier 
transform  are  calculated  on  some  circular  arcs.  Therefore,  prior  to  the  inverse 
FFT  operation,  an  interpolation  technique  should  be  used  to  recover  the  values 
of  the  function  over  a  square  grid.  An  interpolation-free  version  of  the 
Synthetic  Aperture  technique  is  proposed  in  [Nahamoo82],  and  a  very  simple 
implementation  is  presented  in  [Nahamoo84].  However,  the  proposed  procedure 
requires  no  interpolation  only  if  the  sampling  interval  on  the  transmit,  the 
receive  line  and  in  the  object  are  all  the  same.  This  constraint  may  be  very 
restrictive  in  a  practical  implementation  of  such  a  system. 


Figure  38 

For  a  fixed  orientation  of  the  object  in  the  transmission  mode,  the  object’s 
Fourier  transform  is  obtained  in  the  region  shown. 
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Figure  3.9 

For  two  rotational  positions  of  the  object,  90  degrees  apart,  the  object's 
Fourier  transform  is  recovered  inside  a  circle  of  radius  v2 k0.  The  cover¬ 
age  in  the  region  C  is  double. 
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3.3.2  Reflection  Mode 


In  this  section,  we  show  that  in  the  Synthetic  Aperture  approach,  instead 
of  rotation  of  the  object,  the  backscattered  field  together  with  the  forward 
scattered  field  can  be  used  to  recover  the  object’s  Fourier  transform  inside  a 
circle  of  radius  2k0,  rather  than  v/2k0  achieved  in  the  transmission  mode  with 
two  views.  Thus  by  having  two  fixed  parallel  transducer  arrays  on  both  sides 
of  the  object  and  no  mechanical  movements,  the  object  cross  section  can  be 
reconstructed.  In  this  approach,  in  addition  to  the  data  collected  in  the 
transmission  mode,  for  all  the  positions  of  the  transmitter  on  both  of  the 
transducer  lines,  the  scattered  field  on  all  positions  of  the  receiver  on  the  same 
line  is  measured.  This  new  set  of  data  provides  a  coverage  in  the  Fourier 
domain  of  the  object,  inside  a  circle  of  radius  2k0,  which  is  the  complement  of 
the  coverage  achieved  in  the  transmission  mode  for  that  view  of  the  object. 

We  propose  two  ways  of  implementing  a  system  working  in  reflection 
mode.  One  approach  is  to  use  two  transducers  for  transmission  and  reception. 
Obviously  the  scattered  field  at  those  points  where  the  transducers  overlap 
physically  can  not  be  measured.  Perhaps  the  field  at  those  points  can  be 
estimated  bv  the  average  of  the  field  at  the  neighboring  points  without  causing 
any  major  distortion  in  the  reconstruction.  A  second  approach  is  to  replace  the 
two  separate  transducers  moving  on  a  line  by  an  array  of  transducers;  however, 
the  cross  coupling  between  the  elements  of  the  array  makes  the  CW 
transmission  too  erroneous.  To  avoid  the  cross-coupling  problem,  one  can 
transmit  a  long  burst  of  CW  pulse,  and  range  gate  the  received  waveform.  The 

2d 

length  of  the  burst,  T,  has  to  be  shorter  than  — ,  where  d  is  the  distance  from 

co 

the  object  to  the  array  as  shown  in  Fig.  3.10,  and  c0  is  the  speed  of 
propagation.  Before  the  first  echo  reaches  the  transducers,  the  transmitter  is 
switched  to  the  receiving  circuitry  and  can  detect  the  scattered  field.  On  the 
other  hand,  since  the  process  should  be  treated  as  a  broadband  process,  the 
measured  scattered  field  has  to  be  the  summation  of  the  scattered  fields  from 


all  the  points  of  (he  object,  and  not  only  from  a  certain  range  of  the  object. 
Therefore,  T,  the  length  of  th  burst, 

2S 

has  to  be  chosen  so  that  it  is  longer  than  — ^  +  t^,  where  S  is  the  size  of  the 

object  (see  Fig.  3.10)  and  tyy  is  the  time  required  for  the  quadrature  receiver  to 
estimate  the  magnitude  and  phase  of  the  signal.  Thus  we  have 
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—  +  tw  <  T  <  — -  (3.48) 
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In  this  expression  d  was  replaced  by  — - — ,  where  D  (=xr-xt)  is  the  distance 

<w 

between  the  arrays.  This  imposes  the  following  restriction  on  the  object  size 

s  <  J  -  J  tw  ,  (3.40) 

which  means  that  the  object  size  should  be  at  least  three  times  smaller  than 
the  distance  between  the  transducer  arrays. 

The  mathematical  steps  in  the  derivation  of  a  relationship  between  the 
reflected  data  and  the  object  function  is  quite  similar  to  that  of  the 
transmission  case.  The  only  difference  appears  in  the  expansion  of  the  Green’s 
function.  Instead  of  expanding  the  Green’s  function  in  the  forward  direction, 
we  expand  it  in  the  reverse  direction 

g(r)T0)  =  -j-  f  -  d/?  ,  X<XQ  ;  (3.50) 

4r  4  a 

a  is  now  defined  as 

a  =  s/kj  -  f  (3.51) 

The  processing  of  the  backscattered  field  measured  on  the  transducer  line 
at  x=xt  results  in  the  coverage  of  the  object’s  Fourier  transform  in  the  region 
shown  in  Fig.  3.11.  We  will  prove  that  this  mapping  is  two-to-one,  or  in  other 
words  the  value  of  the  object's  Fourier  transform  at  any  point  in  the  region 
shown  in  Fig.  3.11  can  be  calculated  from  two  and  only  two  points  of  the 
processed  data  u3.  When  the  other  transducer  array  at  x=xr  is  used,  the 
coverage  area  would  be  the  symmetric  image  of  the  area  shown  in  Fig.  3.11 
with  respect  to  the  ky  axis.  Therefore,  the  object’s  Fourier  transform  would  be 
recovered  inside  a  circle  of  radius  2k0,  when  both  of  the  transducer  lines  are 
used  in  relflection  mode,  in  addition  to  the  normal  transmission  mode  for  that 
view  of  the  object.  In  Fig.  3.12  we  have  shown  this  coverage.  Notice  that  the 
mapping  between  the  data  and  the  object’s  Fourier  transform  is  one-to-one  in 
the  transmission  mode.  However,  repeating  the  data  collection  procedure  for 
the  transmission  mode  by  switching  the  role  of  the  transmit  and  the  receive 
lines  would  make  it  a  two-to-one  relationship  as  well. 
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Figure  3.11 

When  the  transmitter  array  in  figure  10.6  is  utilized  only  in  the  reflection 
mode,  for  a  fixed  orientation  of  the  object,  the  object's  Fourier  transform 
is  recovered  inside  the  region  shown. 
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In  the  following,  we  show  that  by  employing  Eqs.  (3.47)  and  (3.51)  in  a 
lossless  medium,  we  can  recover  the  Fourier  transform  of  the  object  function 
inside  the  circle  of  radius  2k0  centered  at  the  origin  and  outside  the  circles  of 
radius  k0  centered  at  (ko,0)  and  (— ko,0),  as  shown  in  figure  3.11.  We  bring  the 
proof  for  one  set  of  the  backscattered  data  (for  the  region  shown  in  Fig.  3.10) 
while  the  derivation  for  the  other  region  follows  from  the  symmetry  of  the 


problem. 

Let  W  =  (u,v)  be  defined  by 

W  =  A  -  R  (3.52) 

Eq.  (3.51)  can  be  written  in  the  following  form  when  kco  is  real. 

a2  +  =  k02  (3.53) 

Ignoring  evanescent  waves,  the  limits  on  a  are  found  by  considering  Eq.  (3.51). 

-  k0  <  a  <  0  (3.54) 

Considering  Eq.  (3.33),  one  can  write  the  following  equation  and  inequality 

kx2  +  ky2  =  k02  ,  (3.55) 

0  <  kx  <  k0  .  (3.56) 

Substituting  Eq.  (3.52)  in  (3.53)  and  (3.54),  we  get 

(u  +  kx)2  +  (v  +  ky)2  =  k02  ,  (3.57) 

and 

-  k0  <  u  +  kx  <  0  .  (3.58) 


Eqs.  (3.57)  and  (3.58)  represent  a  semicircle  of  radius  kD  centred  at  point 
(— kx ,— ky).  We  shall  denote  this  semicircle  by  SC(  -kx  ,  -ky  ).  On  the  other 
hand  Eqs.  (3.55)  and  (3.56)  represent  a  semicircle  of  radius  k0  centered  at  the 
origin.  This  semicircle  that  contains  all  the  centers  of  the  SC(-kx,-ky)  is 
depicted  by  the  dashed  arc  in  figure  3.13.  A  typical  SC(-kx,-ky)  is  shown  by 
the  solid  arc  in  this  figure.  The  union  of  SC(-kx,-ky)  for  all  possible  center 
locations  (-kx,-ky)  will  provide  a  coverage  of  F(W)  in  some  region  in  the 
Fourier  domain. 

Consider  the  semicircle  SC(-kx,-ky)  shown  by  the  solid  line  in  Fig.  3.13. 
Let  A(kx,ky)  be  the  center  of  this  semicircle  and  D(u,v)  a  point  on  this 
semicircle.  From  equation  (3.58)  we  have  u  <~kx.  Since  kx  >  0,  u  has  to  be 
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Figure  3.13 

In  the  reflection  mode,  the  object’s  Fourier  transform  is  obtained  over  the 
semi-circles  SC(-kt-kf),  one  of  which  is  the  semi-circle  centered  at  point 
A .  By  moving  this  point  on  the  dashed  serai-circle,  the  coverage  shown  in 
figure  10.11  is  achieved. 
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negative,  thereby  the  coverage  is  restricted  to  the  region  u  <  0.  In  the  triangle 
OAD,  we  have 

OD  <  AD  +  OA  ,  (3.59) 

or  in  terms  of  u  and  v 

u2  +  v2  <  (2k0)2  .  (3.60) 


This  inequality  shows  that  the  coverage  will  be  inside  the  circle  centered  at.  the 
origin  with  the  radius  of  2k0. 

Let  the  point  B  be  the  edge  of  the  semicirlce  SC(-kx,-ky)  in  the  region 
v  <  0.  First  we  consider  the  third  quadrant  (  v<0  )  and  make  use  of  the 
symmetry  of  the  problem  with  respect  to  the  u  axis  and  extend  it  to  the  second 
quadrant  (  v>0  )  as  well.  It  is  clear  from  Fig.  3.13  that  the  angle  DOE  is 
always  less  than  or  equal  to  the  angle  BOE.  The  relation  between  the  tangent 
of  these  two  angles  can  be  written  as 

tan(DOE)  >  tan(BOE)  .  (3-61) 

The  angle  BOE  is  given  by 


BOE  = 


APE 

2 


(3.62) 


We  also  know  that 


tan!  AQE  1  =  t  ~  cos(AOE) 

1  2  1  sin(AOE) 

Using  equations  (3.61),  (3.62)  and  (3.63),  we  obtain 
sin(AOE)  tan(DOE)  >  1  —  cos(AOE)  . 


(3.63) 

(3.64) 


Replacing  AOE  and  DOE,  by  their  expressions  in  terms  of  kx,  ky,  u  and  v,  we 
get 


K 

K 


_u 

V 


(3.65) 


Making  use  of  Eq.  (3.57),  the  inequality  above  is  rewritten  as 

u2  +  (v  +  k0)2  >  k02  (3.66) 

The  inequalities  (3.60)  and  (3.66)  define  the  Fourier  domain  coverage  in  the 
region  u  <  0  ,  v  <  0  to  be  outside  the  circle  centered  at  (0,-ko)  with  radius  k0 
and  inside  the  circle  at  the  origin  with  radius  k0. 


Since  the  problem  is  symmetric  with  respect  to  the  u  axis,  one  can  easily 
generalize  the  result  achieved  in  the  third  quadrant  to  the  second  quadrant.  So 
the  inequality  (3.66)  for  the  whole  half  plane  u  <  0  is  written  as 

u2  +  (|  v|  +  k0)2  >  k02  u  <  0  (3.67) 

Finally  the  coverage  area  is  defined  by  the  inequalities  of  equation  (3.60)  and 
(3.67)  and  is  depicted  by  the  shaded  area  shown  in  figure  3.11. 

Up  to  this  point,  we  have  proved  that  any  data  point  maps  into  the  region 
shown  in  figure  3.11.  To  prove  that  the  coverage  is  indeed  that  region,  we  also 
have  to  prove  that  any  point  in  that  region  belongs  to  a  semicircle  SC(-kx,-ky) 
for  some  kx  and  ky. 

In  the  following,  we  prove  that  indeed  for  any  point  in  the  shaded  area 
shown  in  figure  3.11,  there  exists  two  and  only  two  semicircles  SC(-kx,-ky). 
We  present  the  proof  only  for  the  third  quadrant.  One  can  make  use  of  the 
symmetry  of  the  problem  with  respect  to  the  u  axis  to  extend  it  to  the  whole 
half  plane  u  <  0.  Let’s  take  the  point  D(u,v)  in  the  coverage  area  where 
v  <  0,  as  shown  in  Fig.  3.14.  The  centers  of  any  circle  of  radius  k0  passing 
through  the  point  D,  has  to  be  located  at  the  intersection  of  a  circle  of  radius 
k0  centered  at  D  and  the  circle  of  radius  k0  centered  at  the  origin,  which  we 
call  C(0,0).  Since  the  distance  between  the  centers  of  these  two  circles  is  less 
than  2k0,  there  always  exists  two  intersection  points.  We  have  to  prove  that 
these  intersection  points  are  always  in  the  region  u  <  0.  and  also  show  that 
the  u  component  of  D  is  smaller  than  that  of  the  two  intersection  points,  so 
that  D  lies  on  the  left  half  of  the  circles  passing  through  it. 

To  prove  that  the  two  intersection  points  A  and  B  lie  in  the  left  half 
plane,  we  show  that  the  distance  of  the  point  D  to  any  point  P  on  the  right 
half  of  C(0,0)  is  larger  than  k0.  The  coordinates  of  point  D(u,v)  and  P(up,vp) 


satisfy  the  following  equations 

u2  +  (v  +  ky)2  >  k02  ,u  <  0  ,  v  <0  and  (3.68) 

up2  +  vP2  =  k02  ,uP  >  0  (3.6ft) 

The  distance  DP  is  given  by 

DP2  =  (u-up)2  +  ( v  v p )2  (3.70) 

Using  Eq.  (3.69),  we  can  rewrite  Eq.  (3.70)  in  the  following  form 

DP2  =  u2  +  (v  +k0)2  -  2uuP  -  2v(k0  +  vP)  (3.71) 


Considering  the  inequality  in  (3.68)  and  the  fact  that  the  last  two  terms  in  the 


Figure  3.14 

The  points  A  and  B  are  the  center  of  semi-circles  SC{-kt-kf)  which 
pass  over  the  point  D.  Point  D  is  restricted  to  the  region  shown. 


equation  above  are  always  positive,  we  obtain 


DP  >  kr 


(3.72) 


Therefore,  the  points  A  and  B  can  not  be  on  the  right  half  of  C(0,0)  and  are  in 
the  half-space  u<0. 

Let’s  call  the  u  components  of  the  intersection  points  A  and  B,  uA  and  uB 
respectively.  It  is  left  to  prove  that  u  is  smaller  than  uA  and  uB.  uA  and  uB 
can  be  expressed  as 

uA  =  -  k0  sin(DOE  -  DOA)  and  (3.73a) 

uB  =  -  k0  sin(DOE  +  DOA)  (3.73b) 

Since  the  point  D  lies  in  the  third  quadrant,  the  angle  DOE  is  limited  to  the 
following  range 


-  >  DOE  >  0  . 

2  - 

We  also  have  the  following  restriction  on  the  angle  DOA 
DOE  >  DOA  >  0 

Using  the  inequalities  in  (3.74)  and  (3.75),  we  get 


—  >  DOE  -  DOA  >  0 

2  ~ 


(3.74) 


(3.75) 


(3.76a) 


2~  >  DOE  +  DOA  >  0  (3.76b) 

Equations  (3.73a,b)  can  be  combined  to  obtain 

uB  =  uA  -  2  k0  cos(DOE)  sin(DOA)  (3-77) 

Considering  the  ranges  of  DOE  and  DOA  (the  inequalities  in  (3.74)  and  (3.75)), 
the  second  term  in  equation  (3.77)  is  always  negative,  thus  we  have  uB  <  uA. 
Thereby  one  only  needs  to  prove  that  u  is  smaller  than  uB.  Equation  (3.73b) 
can  be  rewritten  as 

nB  =  -  k„  (  sin(DOE)  cos(DOA)  -  cob(DOE)  sin(DOA)  ) 


(3.77) 


-  k„  Mn(l)OE)  cos(DOA)  1  + 


tan(DOA) 


tan(DOE) 


-  u  OD 


j  +  tan(DOA) 
tan(DOE) 


=  -  k„  ■= 


0  OD  2k0 


tan(DOA) 

tanjDOE) 


From  the  inequalities  in  Eqs.  (3.74)  and  (3.75),  we  have 
tan(DOE)  >  tan(DOA)  >  0 
This  inequality  can  be  written  into  the  following  form 


1  >  1 
“  2 


1  + 


tan(DOA) 


tan(DOE) 


>  1 
~  2 


(3.78) 


(3.79) 


(3.80) 


The  two  inequalities  in  (3.78)  and  (3.80)  show  that  u  <  uB.  The  proof  is 
complete. 


3.3.3  Lossy  Medium 

The  theory  presented  in  section  3.3.1  is  general,  in  the  sense  that  it 
includes  the  case  of  an  attenuating  as  well  as  lossless  medium.  Consider  an 
attenuating  medium  with  the  complex  wavenumber  kc0.  The  parameters  kx 
and  a,  as  given  by  Eqs.  (3.33)  and  (3.40),  become  complex,  and  as  a  result  the 
first  argument  of  F  in  Eq.  (3.47)  (a~kx)  also  becomes  complex.  This  means 
that  by  processing  the  data  according  to  Eq.  (3.47),  one  can  find  the  values  of 
the  object  “Fourier-Laplace”  transform  over  a  certain  region.  By  “Fourier- 
Laplace”  transform  of  a  two-dimensional  function  we  mean  Laplace  transform 
in  one  direction  and  Fourier  transform  in  the  perpendicular  direction  (see  Eq. 
(3.30)). 

Let’s  take  the  vectors  p  and  t[  in  Eq.  (3.30)  to  be  represented  by 

p  =  ur  +  v  ,  (3.8ta) 

cj  =  U;  ,  (3.81b) 

where  u  =  ur  +  j  U;  is  now  a  complex  variable  (see  Eq.  (3.53)).  Let's  define 
W  =  (u.v)  as 

W  =  A  -  Fv  (3.82) 


Eq.  (3.40)  can  be  rewritten  as 


O 

a“ 


+  r- 


(3.83) 


Substituting  Eq  (3.82)  in  (3.83),  we  have 


(u  +  kj-  +  (v+ky)-  -  kCg  .  (3.84) 

On  the  other  hand,  kx  and  ky  are  related  by  Eq.  (3.33). 

k,2  +  k,2  =  k„2  (3.85) 

Eqs.  (3.84)  and  (3.85)  can  be  solved  for  u;.  The  derivation  can  be  found  in 
appendix  D,  the  result  is 

76Uj5  +  05ui4  +  ^4ui3  +  73Uj2  +  72Ui  +  7i  =  0  ,  (3.86) 


where 


7s  ~  -ur 


75  -  4k0e0 


74  =  "2ur(ur2  -  v2) 


73  =  -4k0€0(‘2ur2  ~  v-)  , 

72  =  -urK4  +  v4  -  4(k02-£02v2  +  2ur2v2  )  , 


(3.87a) 

(3.87b) 

(3.87c) 

(3.87d) 

(3.87e) 


7i  =  4ur2k0e0(ur2  +  v2)  .  (3.87f) 

The  satisfactory  solution  of  U;  should  result  in  positive  real  values  for  a  and  kx 
for  the  transmission  mode.  Eq.  (3.86)  defines  the  coverage  surface  in  the 
three-dimensional  space  (ur,uitv). 

Due  to  the  complexity  of  Eq.  (3.86),  one  has  to  take  recourse  to  numerical 
approaches  to  find  the  solution) s)  of  U;.  To  prove  or  disprove  that  the  mapping 
to  the  three-dimensional  space  (ur,Uj,v)  is  one-to-one  seems  to  be  very  involved 
if  not  impossible.  By  applying  the  well  known  Newton's  method,  we  solved  Eq. 
(3.85)  for  values  of  ur  and  v  over  a  square  grid  numerically  and  found  the 
mapping  to  be  a  one-to-two  relationship.  However,  because  of  the  round  off 
errors  and  the  small  number  of  points  (over  the  square  grid)  examined  we  can 
not  generalize  our  conclusion  from  the  numerical  results.  Choosing  the  solution 
with  the  smallest  magnitude,  we  have  plotted  the  corresponding  coverage 
surface  in  Fig.  3.15,  for  a  typical  attenuation  constant  of  3.0  db/cm.  Those 
values  of  u;  with  their  magnitude  larger  than  the  maximum  scale  in  Fig.  3.15 
were  set  to  zero.  As  the  attenuation  goes  to  zero,  the  surface  in  Fig.  3.15  will 
collapse  on  the  urv  plane  and  the  coverage  surface  will  be  the  same  as  that 
shown  in  Fig.  3.8,  as  expected.  The  coverage  surface  in  Fig.  3.15  makes 
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Figure  3.15 

The  curved  surface  shows  the  coverage  achieved  in  the  object’s  “Fourier* 
Laplace”  domain  when  the  medium  is  attenuating.  The  attenuation  of  the 
medium  was  taken  to  be  3.0  db/cm. 


impossible  the  inversion  process  to  the  object  domain.  To  apply  the  inverse 
Laplace  transform  techniques  the  coverage  surface  has  to  be  parallel  to  the  urv 
plane  [Krylov69,  Krylov77,  Beilman66j. 

A  possible  way  of  avoiding  the  “Fourier-Laplace”  transform  operation  is 
to  impose  the  following  constraint  on  the  processed  data  [Nahamoo82] 

|*|=|kx|  .  (3.88) 

Then  the  coverage  area  would  be  limited  to  the  v  axis  (  |  vj  <  2k0  ),  instead  of 
the  region  shown  in  Fig.  3.8.  By  rotating  the  object  180  degrees,  the  object’s 
Fourier  transform  can  be  recovered  over  radial  lines  (see  Fig.  3.16)  and  the 
object  can  be  reconstructed  by  means  of  Fourier  inversion.  As  a  matter  of 
fact,  the  limiting  constraint  defines  the  intersection  of  the  coverage  surface 
shown  in  Fig.  3.15  with  the  plane  urv.  This  is  the  only  reconstruction 
technique  available  today  for  an  attenuating  medium;  however,  the  amount  of 
data  required  makes  it  impractical  with  the  present  technology. 
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CHAPTER  4 

A  SIMULATION  STUDY  OF 
MULTIPLE  SCATTERING  AND  ATTENUATION 
IN  DIFFRACTION  TOMOGRAPHY 


All  Diffraction  Tomography  algorithms  available  today,  assume  the  first 
order  Born  or  Rytov  approximations;  which  implies  that  they  are  only  good  for 
reconstructing  weakly  scattering  objects.  In  a  multi-component  object  this 
translates  into  ignoring  the  interaction  between  the  components.  However  in 
most  cases  of  practical  interest  the  object  would  be  composed  of  several 
distinct  components.  Our  objective  in  section  4.1  is  to  study  the  effect  of  the 
interaction  between  object  components  on  the  reconstructed  object,  in  both  the 
conventional  and  the  Synthetic  Aperture  techniques.  In  this  section,  we  have 
also  looked  into  different  approaches  to  reduce  the  multiple  scattering  effects  in 
the  conventional  technique. 

Although  the  attenuation  in  tissue  is  not  severe  for  ultrasound  fields,  for 
microwave  it  is  on  the  order  of  3  db/cm  at  4  GHz.  Since  there  exists  no 
practical  reconstruction  technique  for  a  lossy  medium  at  the  present  stage,  it  is 
of  importance  to  study  the  performance  of  the  current  reconstruction 
techniques  in  the  presence  of  attenuation.  We  have  presented  this  study  in 
section  4.2. 

The  data  (projections)  required  for  the  computer  simulations  in  this 
chapter  were  provided  by  the  algorithms  discussed  in  chapter  2. 

4.1  Artifacts  Caused  by  Multiple  Scattering  Phenomena 

As  discussed  in  chapter  9,  the  scattering  interaction  between  the 
components  of  an  object  can  cause  the  scattered  field  on  the  receiver  line  to 
become  quite  different  from  the  summation  of  the  first  order  scattered  fields. 
First,  order  fields  are  calculated  by  adding  the  contributions  from  each 
component  while  the  others  are  assumed  to  be  absent.  Since  the  reconstruction 
algorithms  are  based  on  the  weak  scattering  assumption,  the  effect  of  higher 
order  scattering  contributions  present  in  the  projection  should  appear  as  a 
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distortion  term  in  the  final  reconstruction.  We  consider  cases  in  which  the 
reconstruction  of  any  of  the  object’s  components  by  its  own  is  satisfactory. 

Both  the  conventional  and  Synthetic  Aperture  techniques  will  be  tested  to 
observe  the  sensitivity  of  each  to  multiple  scattering  phenomena.  In  section 
4.1.2,  we  consider  several  possible  approaches  to  reduce  the  effect  in  the 
conventional  technique. 

4.1.1  Artifacts  in  Conventional  Diffraction  Tomography 

The  computing  procedure  discussed  in  chapter  9  was  used  to  generate  a 
certain  number  of  projections  over  360  degrees  around  the  object.  Each 
projection  is  the  scattered  field  measured  on  a  receiver  line  of  finite  size  when  a 
multi-component  object  is  illuminated  by  a  plane  wave.  It  should  be  noted 
that  there  has  been  a  huge  computational  effort  required  to  generate  the 
projections.  For  example,  generating  64  projections  of  an  object  consisting  of 
three  cylinders  took  more  than  a  couple  of  hours  of  CPU  time  on  the  FPS- 
AP120B  array  processor. 

In  all  the  simulations  presented  in  this  section,  64  projections  were 
calculated  around  the  object  and  a  bilinear  interpolation  based  algorithm  was- 
then  used  to  reconstruct  the  object  cross  section;  the  result  follows. 

In  the  reconstructions  shown  in  Figs.  4. 1-4.4,  we  have  shown  the 
magnitude  of  the  deviation  of  the  reconstructed  refractive  index  from  that  of 
the  background,  which  was  assumed  to  be  unity.  Plots  labeled  (a)  show  the 
reconstruction  obtained  when  the  projections  were  generated  by  ignoring  the 
second  and  higher  order  scattered  fields.  On  the  other  hand,  the  plots  labeled 
(b)  show  the  reconstructions  when  the  second  order  fields  are  included  in  the 
projections. 

In  Fig.  4.1  the  change  in  the  refractive  index  of  the  6\  cylinders  was  set  to 
2  percent.  Although  there  is  some  distortion  introduced  in  the  direction  of  the 
line  joining  the  center  of  the  cylinders,  it  is  negligible.  However,  when  the 
refractive  index  change  is  increased  in  Fig.  4.2  to  3  percent,  the  distortion 
becomes  quite  noticeable;  and  in  Fig.  4.3  a  5  percent  change  in  refractive  index 
is  enough  to  cause  the  distortion  to  dominate  the  reconstruction.  When  the 
number  of  cylinders  is  increased,  the  distortion  is  higher,  as  seen  in  Fig.  4.4. 
This  is  expected  because  in  this  case  there  are  more  projections  affected  by 
second  order  scattering. 
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Figure  4.1 

Cross  section  of  a  two-component  object  reconstructed  using  the  conven¬ 
tional  diffraction  reconstruction  technique.  Diameter  of  the  cylinders  is  6 
X  and  their  refractive  index  is  1.02.  a)  Only  the  first-order  scattered  fields 
are  used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scat¬ 
tered  fields  are  included  in  the  projection  data  for  this  reconstruction. 
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Figure  -4. 2 

Cross  section  of  a  two-component  object  reconstructed  using  the  conven¬ 
tional  diffraction  reconstruction  technique.  Diameter  of  the  cylinders  is  6 
\  and  their  refractive  index  is  1  03.  a)  Only  the  first-order  scattered  fields 
are  used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scat¬ 
tered  fields  are  included  in  the  projection  data  for  this  reconstruction. 


Figure  4.3 

Cross  section  of  a  two-component  object  reconstructed  using  the  onven- 
tional  diffraction  reconstruction  technique.  Diameter  of  the  cylinders  is  6 
X  and  their  refractive  index  is  1.05.  a)  Only  the  first-order  scattered  fields 
are  used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scat¬ 
tered  fields  are  included  in  the  projection  data  for  this  reconstruction. 
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Figure  41 

Cross  section  of  a  three-component  object  reconstructed  using  the  conven¬ 
tional  diffraction  reconstruction  technique.  Diameter  of  the  cylinders  is  6 
X  and  their  refractive  index  is  1.05.  a)  Only  the  first-order  scattered  fields 
are  used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scat¬ 
tered  fields  are  included  in  the  projection  data  for  this  reconstruction. 
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4.1.2  Different  Approaches  to  Reduce  the  Artifacts 

We  have  experimented  with  three  different  approaches  for  reducing  the 
distortions  caused  by  multiple  scattering  in  the  conventional  Diffraction 
Tomography  technique.  First  a  statement  about  the  motivation  for  this 
experiment:  in  chapter  9,  we  showed  that  even  when  the  object  components  are 
weakly  scattering,  the  multiple  scattering  effects  become  important  if  the 
components  are  blocking  each  other.  By  blocking  we  mean  the  line  joining  the 
component  centers  is  parallel  with  the  propagation  vector  (see  Fig.  4.5a).  In 
tomographic  reconstructions,  these  multiple  scattering  effects  exhibited 
themselves  as  distortions,  especially  along  the  lines  joining  the  components.  In 
the  first  approach  for  reducing  these  distortions,  we  use  oblique  receptions  for 
all  projections.  In  this  method,  the  receive  array  is  made  non-parallel  with  the 
transmit  array,  as  shown  in  Fig.  4.5b.  In  this  way  the  blocking  effect,  which  is 
quite  severe  when  the  propagation  vector  and  the  line  joining  the  object 
components  (for  a  two-component  object)  are  collinear,  is  eliminated.  /Vs  will 
be  demonstrated  in  section  4.2. 1.1,  although  this  arrangement  substantially 
reduces  the  distortions  in  spaces  between  the  object  components,  the  price  paid 
is  a  small  increase  in  the  general  distortion  all  through  the  object.  This 
approach  seems  appropriate  for  two-component  objects. 

In  a  second  approach  intended  for  more  complex  objects,  the  object  was 
illuminated  with  several  non-collinear  plane  waves  for  each  projection,  for  each 
rotational  position  of  the  receive  array  (see  Fig.  4.6).  Each  plane  wave 
generates  a  complete  frequency  domain  coverage;  these  are  averaged  prior  to 
Fourier  inversion  for  object  reconstruction.  The  rationale  for  this  method  is 
that  the  errors  caused  by  blocking  vary  for  the  three  different  illuminating 
waves,  one  reason  being  that  some  of  them  correspond  to  the  oblique  reception 
concept  mentioned  above.  Therefore,  averaging  in  the  frequency  domain 
should  diminish  the  effect  of  the  errors  caused  by  multiple  scattering.  To 
substantiate  this,  we  have  shown,  in  section  4.1.2.!,  reconstructions  obtained 
using  this  idea. 

We  have  also  considered  another  approach  which  is  similar  in  spirit  to  the 
.Synthetic  Aperture  technique.  For  each  position  of  the  receiver  line  the  object 
is  illuminated  by  incident  plane  waves  from  several  direct  ion-1.  (However,  unlike 
the  second  approach,  no  averaging  is  performed  in  the  frequency  domain.)  The 
object  is  now  rotated  by  a  large  angle  to  \  e  ld  a  similar  coverage  over  areas  left 
blank  before:  this  procedure  is  repeated  until  no  holes  are  left  in  the  frequence 
domain.  The  rationale  here  is  that  with  each  plane  wave  at  a  different 
orientation  with  respect  to  the  receiver  line,  the  severity  of  the  blocking  effect 
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Figure  4.5 

(a)  The  condition  when  the  object  components  block  each  other  is  illus¬ 
trated  here. 

(b)  With  the  oblique  reception  illustrated  here  the  "direct  blockage"  is 
avoided.  Direct  blockage  occurs  when  the  direction  of  the  propagation  of 
the  incident  field  is  lined  up  with  the  line  joining  the  two  object  com¬ 
ponents  and  the  normal  to  the  receiver  array. 


Figure  4.6 

In  the  second  approach  several  plane  waves  are  used  to  illuminate  the  ob¬ 
ject.  The  frequency  domain  coverage  obtained  with  the  three  plane  waves 
shown  here  are  averaged  prior  to  inversion. 
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caused  by  multiple  scattering  should  be  reduced.  In  this  case,  the  amount  of 
data  required  is  close  to  that  for  the  usual  diffraction  tomography.  It  was 
found  that  this  approach  does  result  in  better  reconstruction  of  multi- 
component  objects. 


4.1. 2.1  First  Approach 

Examine  the  reconstruction  of  a  two  component  object  in  Fig.  4.3;  there 
the  blocking  effect  led  to  distortion  that  was  relatively  large  in  regions  between 
the  components.  It  was,  therefore,  felt  that  if  we  instead  used  for  each 
projection  the  oblique  reception  scheme  shown  in  Fig.  4.5b,  that  would 
eliminate  the  blocking  effect,  at  least  for  regions  between  the  components.  For 
the  projection  angle  shown  in  Fig.  4.5b  and  with  this  approach,  the  first  order 
fields,  Ujj  and  u2o,  from  each  component  can  propagate  unperturbed  to  the 
receiver  line,  i.e.  without  any  substantial  interaction  with  the  other  component 
Although  the  direct  blocking  (the  components  being  in  each  other’s  shadow  and 
also  being  lined  up  with  the  direction  of  the  incident  plane  wave)  is  eliminated 
by  this  method,  this  leads  to  reduced  reconstruction  distortion  in  spaces 
between  the  components.  Many  more  projections  experience  multiple 
scattering  now,  albeit  at  a  lower  level.  While  for  the  direct  reception  method 
there  is  only  one  projection  where  the  normals  to  the  transmit  array,  the 
receive  array,  and  the  line  joining  the  centers  of  the  object  components  are  all 
lined  up,  for  the  oblique  reception  case  there  will  be  more  directions  in  which 
the  line  joining  the  object  centers  will  be  perpendicular  to  either  the  receive 
array  or  the  transmit  array,  ,4s  a  consequence,  the  oblique  reception  method 
reduces  the  reconstruction  distortion  between  the  components,  but  it  raises 
somewhat  the  overall  level  of  distortion.  This  is  illustrated  in  the  section  4. 1.2.4 
with  computer  simulations. 

With  regard  to  the  coverage  generated  in  the  Fourier  domain  with  the 
oblique  reception  approach,  it  is  unaffected.  When  the  incident  plane  wave  is 
normal  to  the  receiver  line,  the  Fourier  transform  of  the  projections  gives  the 
values  of  the  Fourier  transform  of  the  object  on  a  semicircle,  as  shown  in  Fig. 
4.7a.  When  the  incident  wave  is  at  an  oblique  angle  to  the  receiver  line,  the 
coverage  is  still  on  a  semi-circle,  but  positioned  differently.  This  can  be  seen  in 
Fig.  1.7b.  Obviously  by  rotating  the  object  from  0  to  360  degrees,  the  arc  0.4 
will  cover  all  the  area  inside  a  circle  of  radius  r,.  where  r(  is 

r,  =  OA  {  1.1 1 
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Figure  4.7 

(a)  When  the  incident  plane  wave  is  normal  to  the  receiver  line,  the 
Fourier  transform  of  the  projection  gives  the  values  of  the  object's  Fourier 
transform  over  a  semi-circle  as  shown  here. 

(b)  When  the  incident  wave  is  at  an  oblique  angle  with  the  receiver  line, 
the  coverage  is  still  on  a  semi-circle,  although  positioned  differently. 
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9  is  the  oblique  angle  of  incidence.  A  smaller  circle  also  will  be  covered  by 
rotation  of  the  smaller  arc  OB;  the  radius  of  this  circle  is  r2  where 

QA°_  A 

r2  =  -  k0  sin(— .  (4.2) 

Therefore  we  get  a  double  coverage  inside  the  smaller  circle. 

4. 1.2. 2  Second  Approach 

The  preceding  method  can  be  made  somewhat  independent  of  the  number 
of  object  components  by  using  several  illuminating  plane  waves  for  each 
angular  position  of  the  receive  array,  as  shown  in  Fig.  4.6,  and  then  averaging 
in  the  Fourier  domain  the  coverages  generated  by  each  of  the  plane  waves. 
The  rationale  for  this  approach  is  that  averaging  in  the  frequency  domain 
would  diminish  in  importance  any  distorted  values  that  result  from  blockages. 
With  several  illuminating  waves,  it  is  unlikely  that  the  received  fields  for  all  of 
them  would  suffer  from  blocking  effects  at  the  same  time. 

In  this  approach  the  coverage  in  the  Fourier  domain  would  be  the  union  of 
the  coverage  areas  for  each  incident  plane  wave.  As  discussed  before  in  the  case 
of  a  single  incident  plane  wave,  there  are  two  circles  associated  with  the 
coverage  in  the  Fourier  domain.  The  points  inside  the  smaller  circle  get  double 
coverage,  while  those  in  the  area  between  the  two  circles  are  covered  only  once. 
It  is  now  clear  that  having  several  incident  plane  waves,  there  will  be  many  of 
these  circles,  precisely  twice  the  number  of  existing  plane  waves.  Therefore,  the 
number  of  times  a  point  at  the  center  of  these  concentric  circles  is  covered 
would  be  equal  to  the  number  of  circles  (twice  the  number  of  incident  plane 
waves).  As  the  point  is  moved  away  from  the  center  of  circles,  getting  out  of 
each  circle  the  multiple  coverage  is  reduced  by  one.  Finally,  it  should  be  noted 
that  the  total  coverage  area  (the  area  in  which  all  the  points  are  covered  at 
least  once)  lies  inside  the  largest  of  the  circles,  with  the  radius  of 

OO°  +  max(0;) 

r  -  2  k„  sin  - - - -  (4.3) 

where  max(0j)  is  the  largest  oblique  angle. 
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4. 1.2. 3  Third  Approach 

In  this  approach,  we  again  use  a  number  of  plane  waves  for  illumination 
for  each  rotational  position  of  the  object.  However,  the  coverages  are  not 
averaged  in  the  frequency  domain.  Instead,  we  use  the  procedure  described 
below. 

If  we  use  five  plane  waves  15°  apart,  the  combined  coverage  in  the 
frequency  domain  for  one  rotational  position  of  the  object  is  as  shown  in  Fig. 
4.8.  The  object  is  now  rotated  by  a  large  angle  to  yield  a  similar  coverage  over 
areas  left  blank  by  the  first  rotational  position,  and  so  on.  With  very  few- 
rotational  positions  of  the  object,  the  complete  frequency  coverage  can  be 
obtained  (The  method  is  similar  in  spirit  to  the  Synthetic  Aperture  approach). 

We  felt  that  this  approach  would  also  yield  some  protection  against  the 
multiple  scattering  distortion  due  to  the  fact  that  for  each  rotational  position 
the  severity  of  the  blocking  effect  would  be  somewhat  reduced  owing  to  the 
plane  waves  being  at  different  angles  to  the  receive  line. 

4. 1.2. 4  Simulation  Results  and  Discussion 

For  the  oblique  reception  method,  we  first  modified  our  computer 
programs  for  projection  data  generation  to  handle  any  arbitrary  direction  of 
the  incident  field.  We  used  this  program  to  generate  the  data  required  for  the 
simulations  done  in  this  section.  Further  programming  had  to  be  done  to  code 
the  reconstruction  algorithms  used  for  the  results  shown  here. 

In  Fig.  4.9a-d  we  have  shown  the  reconstruction  of  a  two-component 
object  with  each  component  of  refractive  index  1.05.  The  components  are  t 
wavelengths  in  diameter  and  12  wavelengths  apart  from  each  other.  The 
receiver  line  is  at  a  distance  of  20  wavelengths  from  the  center  of  coordinates. 
The  sampling  interval  on  this  line  is  0.45  wavelengths  and  there  are  64  points 
in  each  projection.  Plots  a.b.c  and  d  in  Fig.  4.9  represent  the  reconstructions 
of  the  object  when  the  direction  of  the  incident  plane  wave  is.  respectively,  0°. 
7.5°,  15°  and  30°  from  the  normal  to  the  receiver  line.  As  we  expected,  the 
larger  the  off  normal  angle,  the  lower  the  distortion  on  the  line  joining  the 
components.  However,  the  general  level  of  background  distortion  in  the 
reconstruction  increases  with  an  increase  in  this  angle,  as  is  evidenced  by  its 
severity  when  the  angle  is  30°.  The  reconstructions  in  plots  1.10b  and  4.10c  for 
7.5°  and  15°  respectively,  show  some  improvements  over  the  normal  angle  case 
in  plot  4.9a. 
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Figure  -1.8 

In  the  third  approach,  again  a  number  of  plane  waves  are  used  to  il¬ 
luminate  the  object,  however  no  independent  frequency  domain  coverages 
are  generated  from  each  of  these  plane  waves.  For  each  rotational  posi¬ 
tion  of  the  object,  the  coverage  generated  by  all  the  illuminating  waves  is 
extracted  first,  and  then  the  object  is  rotated  till  the  frequency  domain  is 
filled  up. 


Figure  4.9 

Reconstructions  of  a  two-component  object  using  the  first  approach.  The 
cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and  are  positioned  at 
a  distance  of  12  X  from  each  other.  The  projections  are  sampled  at  64 
points  with  a  sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from 
the  center  of  the  coordinates.  The  angle  between  the  direction  of  illumina¬ 
tion  and  the  receiv'.;  line  is:  (a)  0*;  (b)  7.5®;  (c)  15®;  (d)  30®. 
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The  results  obtained  with  the  second  approach  are  illustrated  in  Fig.  4.10. 
Here  the  frequency  domain  coverages  for  each  of  the  plane  waves  shown  in  Fig. 
4.6  are  averaged  prior  to  Fourier  inversion.  Three  plane  waves  are  considered 
at.  angles  of  -15°,  0°  and  15°  in  Fig.  4.11a,  and  -30°,  -15°,  0°,  15°  and  30°  in 
plot  4.11b.  Some  improvement  over  the  first  approach  is  observed,  especially 
in  plot  4.10a  where  the  edges  are  enhanced. 

The  reconstructions  in  Fig.  4.11  were  generated  using  the  third  approach, 
where  we  fill  a  part  of  the  Fourier  domain  for  each  receiver  position  and  several 
incident  plane  waves,  then  rotate  the  object  and  do  the  same  thing  until  a 
complete  coverage  is  obtained.  In  plot  4.11a,  we  have  considered  three  incident 
plane  waves  at  -7.5°  and  7.5°,  and  32  receiver  positions  were  used.  In  plot 
4.11b  the  number  of  plane  waves  was  increased  to  5  with  illumination  angles 
between  -15°  and  15°.  We  see  a  slight  improvement  in  the  reconstructions, 
especially  in  plot  4.11b.  When  the  range  of  angles  is  increased  to  (-30°,  30°), 
the  distortion  in  the  direction  of  the  line  joining  the  cylinders  is  reduced  but  we 
see  an  increase  in  the  general  background  distortion. 

In  the  second  set  of  simulations,  we  changed  the  distance  between  the 
components  to  8  wavelengths  and  their  refractive  index  to  1.08.  Figs.  4.12a-d 
show  reconstructions  of  the  object,  when  the  incident  plane  wave  directions  are, 
respectively,  0°,  10°,  20n  and  30°  from  the  normal  to  the  receiver  line.  We  see 
the  same  changes  in  the  case  where  the  components  were  farther  from  each 
other.  The  interesting  case  is  seen  in  plot  4  12(1.  in  which  the  distortion  in  the 
region  between  the  components  is  significantly  reduced,  but  of  course  some 
distortion  has  been  introduced  all  over  the  reconstruction  area.  Since  the 
components  are  close  to  each  other  in  this  case,  the  receiver  line  has  to  be  quite 
a  bit  oblique  in  this  case  so  that  the  components  may  not  block  each  other 
when  they  are  in  line  with  the  direction  of  propagation  of  the  illuminating  field. 
In  the  reconstruction  obtained  with  the  second  approach  shown  in  Fig.  4.13  the 
recovered  Fourier  transforms  of  the  object  for  different  angles  of  incidence  are 
averaged.  In  plots  a,  b  and  c  there  are,  respectively,  3  incident  plane  waves  in 
the  range  (-10°  ,  10°),  4  in  the  range  (-20°  .  20°)  and  7  in  the  range  (-30°  , 
30°).  In  plot  4.13c,  the  improvement  is  quite  noticeable.  We  have  shown  the 
results  of  the  third  approach  in  Fig.  4.14.  In  plot  4.14a.  there  were  13  receiver 
positions  and  5  incident  plane  waves  in  range  (-20°  ,  20°)  for  each  of  the 
receiver  positions.  Comparing  these  plots  with  1.13a,  we  see  no  improvement. 
As  soon  as  the  range  of  incident  angles  is  increased  to  (-30°  ,  30°)  tn  plot  4.14b. 
the  distortion  between  the  cylinders  is  reduced.  These  results  lead  us  to  the 
conclusion  that  for  closer  spacing  of  components,  the  angle  of  off-normal 
incident  directions  should  become  larger.  We  also  see  a  little  improvement  in 
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Figure  4.10 

Reconstructions  of  a  two-component  object  using  the  second  approach. 
The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and  positioned 
at  a  distance  of  12  X  from  each  other.  The  projections  are  sampled  at  64 
points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from 
the  center  of  the  coordinates,  (a)  3  incident  plane  waves,  with  the  angle  of 
illumination  between  -15*  and  15*.  (b)  3  incident  plane  waves,  with  the 
angle  of  illumination  between  -30*  and  30*. 


Figure  4.11 

Reconstructions  of  a  two-component  object  using  the  third  approach.  The 
cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and  positioned  at  a 
distance  of  12  X  from  each  other.  The  projections  are  sampled  at  64 
points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from 
the  center  of  the  coordinates,  (a)  3  incident  plane  waves,  with  the  angle  of 
illumination  between  -7.5°  and  7.5®  and  32  receiver  positions,  (b)  5  in¬ 
cident  plane  waves,  with  the  angle  of  illumination  between  -15*  and  15® 
and  32  receiver  positions,  (c)  5  incident  plane  waves,  with  the  angle  of  il¬ 
lumination  between  -30®  and  30*  and  16  receiver  positions. 
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Figure  4.12 

Reconstructions  of  a  two-component  object  using  the  6rst  approach  con¬ 
sisting  of  oblique  reception.  The  cylinders  of  refractive  index  1.08  are  4  X 
in  diameter,  and  positioned  at  a  distance  of  8  X  from  each  other.  The  pro¬ 
jections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X.  The 
receiver  line  is  20  X  from  the  center  of  the  coordinates.  The  angle  between 
the  direction  of  illumination  and  the  receiver  line  is:  (a)  0”;  (b)  10”;  (c) 
20”,  (d)  .50”. 
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Figure  4.13 

Reconstructions  of  a  two-component  object  using  the  second  approach. 
The  cylinders  of  refractive  index  1.08  are  4  X  in  diameter,  and  positioned 
at  a  distance  of  8  X  from  each  other.  The  projections  are  sampled  at  64 
points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from 
the  center  of  the  coordinates,  (a)  3  incident  plane  waves,  with  the  angle  of 
illumination  between  -10°  and  lO”.  (b)  5  incident  plane  waves,  with  the 
angle  of  illumination  between  -20'  and  20'.  (c)  7  incident  plane  waves, 
with  the  angle  of  illumination  between  -30"  and  30<>. 


Figure  4.14 

Reconstructions  of  a  two-component  object  using  the  third  approach.  The 
cylinders  of  refractive  index  1.08  are  4  X  in  diameter,  and  positioned  at  a 
distance  of  8  X  from  each  other.  The  projections  are  sampled  at  64  points 
over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  20  X  from  the 
center  of  the  coordinates,  (a)  5  incident  plane  waves,  with  the  angle  of  il¬ 
lumination  between  -20®  ana  20®  and  16  receiver  positions,  (b)  7  incident 
plane  waves,  with  the  angle  of  illumination  between  -30®  and  30®  and  8 
receiver  positions. 


the  results  from  the  second  approach  over  that  of  the  third,  but  the  amount  of 
data  required  for  the  second  one  is  much  higher. 

We  also  tested  the  different  techniques  for  a  three-component  object;  the 
results  are  shown  in  Fig.  4.15-17.  There  is  just  one  component  added  to  the 
geometry  of  the  first  simulation.  Plots  16a-d  show  the  reconstruction  of  the 
object  when  the  incident  plane  wave  direction  is,  respectively,  0°,  10°,  20°  and 
30°  from  the  normal  to  the  receiver  line.  The  improvement  in  the 
reconstruction  is  not  as  noticeable  as  in  the  case  for  a  two-component  object, 
but  still  the  same  effects  are  seen.  The  same  thing  can  be  said  about  the 
results  of  the  second  approach  shown  in  the  next  figure.  In  the  reconstruction 
shown  in  Fig.  4.16.  the  recovered  Fourier  transform  of  the  object  for  different 
angles  of  incidence  is  averaged.  In  plo  s  a,  b  and  c  there  are,  respectively,  3 
incident  plane  waves  in  the  range  (-10°  .  10°),  5  in  the  range  (-20°  ,  20°)  and  7 
in  the  range  (-30°  ,  30°).  We  have  shown  the  results  of  the  third  approach  in 
Fig  4.17.  In  plot  4.17a,  there  were  16  receiver  positions  and  5  incident  plane 
waves  in  the  range  (-20°  ,  20°)  for  ea<  h  of  the  receiver  positions.  When  the 
range  of  angles  is  increased  to  (-30°  ,  30°)  in  plot  4.17b,  the  distortion  in  the 
regions  between  the  components  is  reduced. 

4.1.3  Artifacts  in  Synthetic  Aperture  Diffraction  Tomography 

For  simulation  purposes,  the  transmitting  transducer  can  be  assumed  to  be 
a  point  source.  Even  if  the  transducer  is  not  an  isotropic  radiator,  the 
deconvolution  process  in  Eq.  (3.47)  (devision  by  At( ky ))  would  make  it  function 
as  a  point  source.  Therefore  by  choosing  a  point  source  as  the  transmitting 
transducer,  one  can  avoid  the  deconvolution  computations  and  the  difficulties 
associated  therewith. 

Similar  to  the  conventional  technique,  the  Synthetic  Aperture  technique 
also  calls  for  the  interpolation  to  a  square  grid  in  the  Fourier  domain.  An 
interpolation-free  version  is  proposed  in  [Nahamoo82],  and  implemented  in 
[NahamooSlJ.  However,  the  proposed  procedure  requires  no  interpolation  only 
if  the  sampling  interval  on  the  transmitter  line,  receiver  line  and  in  the  object 
are  all  the  same.  Taking  recourse  to  the  Fourier  domain  interpolation 
technique,  the  values  of  the  object's  Fourier  transform  are  obtained  only  over  a 
finite  number  of  semicircles  in  the  Fourier  domain,  which  calls  for  some  kind  of 
interpolation  scheme.  The  approach  used  in  [Nahamoo82]  was  to  find  the 
neighboring  available  values  of  ky,  assume  the  same  value  of  3  for  each  of  the 
neighboring  ky  values  and  find  its  neighboring  available  values  of  3.  In  Figs. 
4.18(a,b),  the  points  A  and  .V  .  B  and  TV,  and  C  and  C  have  the  same  values 
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Figure  4.15 

Reconstructions  of  a  three-component  object  using  the  first  approach  con¬ 
sisting  of  oblique  reception.  The  cylinders  of  refractive  index  1.05  are  4  X 
in  diameter,  and  positioned  at  a  distance  of  6  X  from  the  center  of  coordi¬ 
nates.  The  projections  are  sampled  at  64  points  over  a  sampling  interval 
of  0.45  X.  The  receiver  line  is  20  X  from  the  center  of  the  coordinates.  The 
angle  between  the  direction  of  illumination  and  the  receiver  line  is  (a)  0*, 
(b)  10*,  (c)  20*,  (d)  30*. 


Figure  -*.16 

Reconstructions  of  a  three-component  object  using  the  second  approach. 
The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and  positioned 
at  a  distance  of  6  X  from  the  center  of  coordinates.  The  projections  are 
sampled  at  64  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line 
is  20  X  from  the  center  of  the  coordinates,  (a)  3  incident  plane  waves, 
with  the  angle  of  illumination  between  -lO9  and  lO®.  (b)  5  incident  plane 
waves,  with  the  angle  of  illumination  between  -20"  and  20*'.  (c )  7  in¬ 
cident  plane  waves,  with  the  angle  of  illumination  between  -30*  and  SO". 


Figure  4.16,  continued. 
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Figure  4.17 

Reconstructions  of  a  three-component  object  using  the  third  approach. 
The  cylinders  of  refractive  index  1.05  are  4  X  in  diameter,  and  positioned 
at  a  distance  of  6  X  from  the  center  of  coordinates.  The  projections  are 
sampled  at  64  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line 
is  20  X  from  the  center  of  the  coordinates,  (a)  5  incident  plane  waves, 
with  the  angle  of  illumination  between  -20*  and  20®  and  16  receiver  posi¬ 
tions.  (b)  7  incident  plane  waves,  with  the  angle  of  illumination  between 
-30®  and  30®  and  16  receiver  positions. 
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of  3,  respectively.  The  values  of  kY  in  Fig.  4.18a  and  4.18b  were  chosen  to  be 
close  to  k0  and  zero,  respectively.  If  a  bilinear  interpolation  scheme  is  used,  the 
four  neighbors  of  a  point  seem  to  be  the  edges  of  a  semi-parallelogram,  as 
expected.  However,  a  nearest  neighbor  scheme  may  result  in  the  completely 
wrong  nearest  neighbor  of  the  point  of  interest  because  of  the  positions  of 
points  A  and  A',  B  and  B' ,  and  C  and  C1  relative  to  each  other.  In  Fig.  4.19 
we  have  shown  a  modified  interpolation  scheme.  Point  C  is  the  point  at  which 
the  value  of  the  object’s  Fourier  transform  is  sought.  After  finding  the 
neighboring  values  of  ky  (points  O,  and  02),  the  intersections  of  the  circle 
going  through  point  C  and  the  two  semicircles  centered  at  O,  and  02  are  found 
(points  C,  and  C2)  and,  finally,  the  neighboring  points  reached  (points  A,.  B,. 
A2  and  Bo). 

The  old  and  the  new  interpolation  techniques  are  compared  in  Figs.  4.20 
and  4.21  for  both  the  nearest  neighbor  and  bilinear  schemes.  We  h  ave  used  the 
computational  technique  discussed  in  chapter  9  to  generate  the  required  data 
for  the  object  reconstruction.  The  number  of  array  elements  in  both  the 
receiver  and  transmitter  has  been  doubled  in  the  simulation  result  in  Fig.  4.21 
compared  to  that  in  Fig.  4.20.  As  expected  there  is  a  noticeable  improvement 
in  the  modified  nearest  neighbor  scheme  compared  to  the  old  nearest  neighbor 
technique.  But  there  seems  to  be  no  change  when  the  bilinear  interpolation 
technique  is  used.  The  results  of  the  old  and  the  new  techniques  become  more 
similar  when  the  number  of  array  elements  is  increased,  which  is  intuitively 
justifiable. 

In  the  formulation  of  the  Synthetic  Aperture  technique,  the  transmit  and 
the  receive  lines  were  assumed  to  possess  infinite  length.  However,  in  the 
simulation  environment,  the  finite  size  of  the  transducer  arrays  may  cause  some 
distortion  in  the  final  reconstruction.  Let's  first  consider  the  receive  line.  For 
those  transducer  points  on  the  edge  of  the  receiver  line  which  are  at  the  far 
field,  the  change  in  the  field  observed  by  moving  to  a  neighboring  point  is  a 
simple  phase  change  due  only  to  the  propagation  effect  of  the  wave.  If  the 
receive  line  is  limited  in  size  to  a  length  several  times  larger  than  the  object's 
size,  the  efTect  would  appear  as  a  lower  energy  at  the  spatial  frequencies  near 
k0,  which  corresponds  to  the  center  frequency  of  the  incident  field.  This 
happens  because  the  high  frequency  oscillations  of  the  field  at  the  end  of  the 
receiver  line  is  simply  zeroed.  Since  k0  is  the  largest  spatial  frequency 
considered  in  the  reconstruction  technique,  thereby  one  can  associate  the  finite 
size  of  the  receive  line  with  the  following  limitation  on  3  in  Fq.  (3.47). 


Figure  4.18 

A  simple  interpolation  scheme  for  the  Synthetic  Aperture  technique 
Points  .4  and  A'  ,  B  and  B  ,  and  C  and  C  lying  on  two  different  semi 
circles  SC(-kj-k  )  posses  equal  values  of  J,  respectively. 

(a)  small  values  of  kt. 


parallel  lines 


(b) 


Figure  4.18,  continued 


Figure  4. 19 

Points  C,  and  C2  are  the  intersection  points  of  the  semi-circles  Ox  and  02 
with  the  circle  centered  at  the  origin  going  through  point  C .  The  neigh¬ 
boring  points  of  C,  and  C\  are  chosen  for  the  kind  of  interpolation 
desired. 
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Figure  4.20 

Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with  diame¬ 
ter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic  Aperture  tech¬ 
nique.  Transmitter  and  receiver  arrays,  positioned  20  X  apart,  have  32  ele¬ 
ments  with  inter-element  distance  of  0  45  X.  (a)  Simple  nearest  neighbor 
interpolation,  (b)  Modified  nearest  neighbor  interpolation,  (c)  Simple  bil¬ 
inear  interpolation,  (d)  Modified  bilinear  interpolation. 
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Figure  4.20,  continued. 
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Figure  1.21 

Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with  diame¬ 
ter  of  6  X  and  refractive  index  of  1.03,  using  the  Synthetic  Aperture  tech¬ 
nique.  Transmitter  and  receiver  arrays,  positioned  20  X  apart,  have  64  ele¬ 
ments  with  inter-element  distance  of  0.45  X.  (a)  Simple  nearest  neighbor 
interpolation,  (b)  Modified  nearest  neighbor  interpolation,  (c)  Simple  bil¬ 
inear  interpolation,  (d)  Modified  bilinear  interpolation. 
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Figure  4.21,  continued. 


|  3\  <  ak0  .  (4.4) 

where  a  is  a  factor  less  than  unity. 

A  similar  argument  can  be  used  for  the  transmit  array,  except  this  time 
the  received  field  for  the  transmitter  positions  at  the  edge  of  the  transmit  array 
are  multiplied  by  a  phase  factor  when  one  moves  to  the  neighboring 
transmitter  position.  A  limitation  similar  to  that  in  Eq.  (4.4)  would  be  imposed 
on  the  magnitude  of  ky 

|  ky  j  <  a  kQ  ( 1.5) 

The  magnitude  of  the  scattered  field  from  a  cylinder  with  refractive  index 
of  1.03  and  diameter  of  6  wavelengths,  as  a  function  of  the  transmitter  and  the 
receiver  positions,  is  plotted  in  Fig.  4.22.  The  effect  of  the  finite  size  of  the 
arrays  is  seen  as  the  chopping  of  the  signal  at  the  end  of  array  positions.  The 
magnitude  of  the  Fourier  transform  of  the  field  with  respect  to  the  transmitter 
and  the  receiver  positions  is  shown  in  Fig.  4.23.  In  plot  (a)  the  number  of 
array  elements  is  32,  while  in  plot  (b)  the  array  elements  have  been  increased 
to  128  elements.  The  clearest  advantage  in  plot  (b)  relative  to  plot  (a)  is  the 
recovery  of  higher  frequencies  at  the  end  of  the  diagonal  line  (higher  values  of 
ky  and  3).  We  expect  a  constant  value  at  those  points  for  which  ky  =  3, 
because  that  corresponds  to  the  value  of  the  object’s  Fourier  transform  at  the 
origin.  To  correct  this  effect  as  observed  in  Fig.  4.23b,  we  normalize  the  result, 
taking  the  value  at  ky  =  3  =  0  as  the  reference  point.  This  point  was  chosen 
because  it  is  the  lowest  frequency  component  available  and,  therefore,  is  less 
distorted. 

It  can  be  easily  shown  that  the  limits  in  Eqs.  (4.4)  and  (4.5)  cause  the 
coverage  in  Fig.  3.8,  to  shrink  to  the  one  shown  in  Fig.  4.24.  Finally  the 
coverage  achieved  for  two  rotational  views  of  the  object  90°  apart,  would  be 
inside  a  circle  of  radius 

r  =  \/2  k0  (  a  -  'Z  l~  a2  )  ;  a>—  ,  (4.6) 


r  =  y/2  k0  (  a  -  V I  -  a2  )  ;  a  >  ,  (4.6) 

rather  than  \/2k0.  To  avoid  the  loss  of  resolution,  one  can  collect  the  data  for 
several  views  of  the  object.  A  typical  coverage  obtained  for  three  rotational 
views  of  the  object  60°  apart  is  shown  in  Fig.  4.25.  The  coverage  in  this  case 
would  be  inside  a  circle  of  radius 
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Figure  4  22 

agnitude  of  the  scattered  field  from  a  lossless  cylinder,  embedded  in  a 
■isless  medium,  as  a  function  of  transmitter  and  receiver  positions.  Di- 
leter  and  refractive  index  of  the  cylinder  are  6  X  and  1.03,  respectively. 
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Figure  -*23 

Magnitude  of  the  Fourier  transform  of  the  scattered  field  from  a  lossless 
cylinder,  embedded  in  a  lossless  medium,  with  respect  to  transmitter  and 
receiver  position.  Diameter  and  refractive  index  of  the  cylinder  are  6  X 
and  1.03,  respectively,  (a)  32  element  transmitter  and  receiver  arrays,  (b) 
128  element  transmitter  and  receiver  arrays. 
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Figure  T25 

When  the  high  frequency  information  in  the  projections  is  lost,  larger 
number  of  views  can  increase  the  resolution  and  decrease  the  distortion  in 
the  reconstructed  image.  The  coverage  shown  corresponds  to  three  rota¬ 
tional  views  of  the  object  60®  apart. 
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In  Fig.  4. ‘26,  we  have  shown  the  reconstructed  Fourier  transform  of  a 
cylinder  with  refractive  index  of  1  03  and  diameter  of  (>  wavelengths.  The 
number  of  array  elements  in  plots  (a)  and  (b)  were  taken  to  be  32  and  12*. 
respectively.  The  severe  distortion  in  the  diagonal  direction  agree*  with  the 
concept  explained  by  Fig.  1.24.  In  Figs.  1.27  and  1.28  we  have  shown  the  real 
(refractive  index)  and  the  imaginary  (difference  in  attenuation  constant)  parts 
of  the  reconstructed  object  function  for  different  sizes  of  the  array*, 
respectively.  Of  course,  the  reconstructed  attenuation  function  should  liava* 
been  zero.  However,  due  to  the  large  size  of  the  object  and  large  value  of  the 
refractive  index,  the  large  estimation  error  in  the  phase  causes  a  transfer  of 
energy  from  the  reconstructed  refractive  index  ter  the  reconstructed  attenuation 
function.  Since  the  medium  attenuation  is  zero,  the  magnitude  of  the 
reconstructed  object  function  would  have  been  a  better  estimate  of  the  object's 
refractive  index.  We  have  chosen  the  real  and  imaginary  parts  only  for  a 
better  two-dimensional  presentation  of  the  distortion.  The  distortion  along  the 
diagonal  direction  in  the  reconstructed  object  is  again  due  to  the  loss  of  the 
high  frequency  components  of  kv  and  J.  One  can  get  rid  of  the  distortion  by 
increasing  the  size  of  the  arrays,  as  shown  in  Figs.  4.27c  and  4.28c.  The 
simulation  results  show  that  the  loss  of  coverage,  shown  in  Fig.  4.24,  is  indeed 
correct  and  a  larger  number  of  elements  in  the  transmitter  and  receiver  would 
solve  the  problem. 

The  simulation  study  presented  in  section  4.1.2  persuaded  us  that  the 
performance  of  the  Synthetic  Aperture  technique  in  the  presence  of  multi- 
component  interactions  would  be  superior  to  that  of  conventional  diffraction 
tomography.  The  programs  developed,  based  on  the  techniques  discussed  in 
chapter  9,  were  used  to  generate  the  required  projections  from  a  multi- 
component  object.  The  computational  time  for  a  three  cylinder  object  was  on 
the  order  of  14  hours  of  CPU  time  in  the  array  processor  FPS-AP120B.  The 
reconstructions  shown  in  Figs.  4.29-31,  for  the  Synthetic  Aperture  technique 
correspond  to  the  ones  presented  in  Figs.  4.2-4  for  the  conventional  technique. 
The  number  of  array  elements  was  taken  to  be  128,  to  avoid  the  effect  of  the 
finite  size  of  the  arrays.  The  arrays  were  positioned  l’>  wavelengths  away  from 
the  center  of  coordinates  and  the  inter-element  distance  was  set  to  0.4) 
wavelengths.  The  results  look  very  interesting  because  there  is  almost  no 
distortion  along  the  line  joining  the  center  of  cylinders  as  was  the  case  in  the 
reconstructions  achieved  using  the  conventional  technique.  The  reason  is  that 
the  blocking  effect  is  minimized  in  this  approach.  Another  interesting  point 
about  the  reconstructions  shown  in  Figs.  4.29-31  is  that  the  components 
blocked  are  more  distorted  compared  to  the  ones  which  are  closer  to  the 
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Figure  4.26 

Magnitude  of  the  reconstructed  object’s  Fourier  transform  in  a  lossless 
medium.  The  object  is  a  lossless  cylinder  with  diameter  and  refractive  in¬ 
dex  of  6  X  and  1.03,  respectively,  (a)  32  element  transmitter  and  receiver 
arrays,  (b)  128  element  transmitter  and  receiver  arrays. 
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Figure  4  27 

Real  part  of  the  reconstructed  object  function  of  a  cylinder  with  diameter 
of  6  X  aDd  refractive  index  of  1.03,  using  the  Synthetic  Aperture  tech¬ 
nique.  Transmitter  and  receiver  arrays,  with  the  inter-element  distance  of 
0.45  X,  are  positioned  20  X  apart,  (a)  32  element  arrays,  (b)  64  element 
arrays,  (c)  128  element  arravs. 
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Figure  1  28 

Imaginary  part  of  the  reconstructed  object  function  (difference  in  attenua¬ 
tion  constant  of  object  and  medium)  of  a  cylinder  with  diameter  of  6  X 
and  refractive  index  of  1.03,  using  the  Synthetic  Aperture  technique. 
Transmitter  and  receiver  arrays,  with  the  inter-element  distance  of  0.45  X, 
are  positioned  20  X  apart,  (a)  32  element  arrays.  ( b )  64  element  arrays, 
(c)  128  element  arrays. 
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Figure  I  29 

Magnitude  of  the  reconstructed  refractive  index  of  a  two-component  ob¬ 
ject,  using  the  Synthetic  Aperture  technique.  The  components  are 
cylinders  with  diameter  of  6  X  and  refractive  index  of  1.03.  Transmitter 
and  receiver  arrays  positioned  30  X  apart,  have  128  elements  with  the 
inter-element  distance  of  0.45  X.  a)  Only  the  first-order  scattered  fu  Ids  are 
used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scattered 
fields  are  included  in  the  projection  data  for  this  reconstruction. 


r  igure  4  .30 

Magnitude  of  the  reconstructed  refractive  index  of  a  two-component  ob¬ 
ject,  using  the  Synthetic  Aperture  technique.  The  components  are 
cylinders  with  diameter  of  6  X  and  refractive  index  of  1.05.  Transmitter 
and  receiver  arrays  positioned  30  X  apart,  have  128  elements  with  the 
inter-element  distance  of  0.-45  X.  a)  Only  the  first-order  scattered  fields  are 
used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scattered 
fields  are  included  in  the  projection  data  for  this  reconstruction. 


Figure  4.31 

Magnitude  of  the  reconstructed  refractive  index  of  a  three-component  ob¬ 
ject,  using  the  Synthetic  Aperture  technique.  The  components  are 
cylinders  with  diameter  of  6  X  and  refractive  index  of  1.05.  Transmitter 
and  receiver  arrays  positioned  30  X  apart,  have  128  elements  with  the 
inter-element  distance  of  0.45  X.  a)  Only  the  first-order  scattered  fields  are 
used  for  generating  the  data  for  this  reconstruction,  b)  Doubly  scattered 
fields  are  included  in  the  projection  data  for  this  reconstruction. 
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transmitter.  Note  that  the  transmitter  array  has  been  positioned  on  the 
negative  side  of  the  axes  x  and  y  in  those  figures.  In  the  case  of  two- 
component  objects  the  effect  is  quite  dear.  In  Fig.  1.51  (threc-oomponet 
object),  the  cylinder  which  is  farther  away  from  the  transmitter  array  in  both 
the  projections,  is  distorted  more  than  the  other  two,  which  seem  to  hold  the 
same  amount  of  distortion. 

4.2  A  Simulation  Study  of  Attenuation  in  Diffraction  Tomography 

Due  to  the  lack  of  any  practical  Diffraction  Tomography  technique  for  a 
lossy  medium,  a  study  of  the  performance  of  the  present  techniques  in  the 
presence  of  attenuation  was  necessary.  Although  in  ultrasonic  imaging  the  role 
of  attenuation  is  minor,  in  microwave  imaging  it  can  not  be  ignored,  because 
microwaves  at  4  GHz  undergo  almost  3  db  of  attenuation  per  centimeter  of 
travel  in  water. 

Let's  consider  the  case  where  the  object  and  the  medium  have  different 
attenuation  constants.  Therefore,  the  wavenumber  in  the  object  (kcfr))  and  the 
medium  (kco)  are  both  complex.  kc(7)  is  defined  as 

kc(T)  =  k(T)  +  j  c(r)  .  (4.8) 

The  reconstruction  algorithms  based  on  the  Fourier  Diffraction  theorem 
estimate  the  complex  object  function  defined  in  Eq.  (3.12).  Since  n(7),  as 
defined  in  Eq.  (3.10)  can  also  be  expressed  in  terms  of  the  wavenumbers  kc(T) 
and  kro  as 


nfr)  -  — —  ,  (4-9) 

kco 

one  can  rewrite  Eq.  (3.12)  in  the  following  form 

Ffr)  =  kc2  -  kco2  .  (4.10) 

Since  kr0  is  known,  both  the  refractive  index  and  the  attenuation  of  the  object 
can  be  estimated  simultaneously. 

Substituting  the  expressions  given  by  Eqs.  (2.1)  and  (4.8)  in  Eq.  (4.10)  and 
restricting  T  to  the  interior  of  the  object,  we  get 

bfr)  =  k02  -  kJ(7)  +  c-(T)  -  cQ2  +  j  2  (  k0  c0  -  kfr)  ((f)  J 

—  -  kn  (  k„  ~  k(T)  j  +  j  2  (  k„  c„  -  k(T)  c(T)  )  (4.11) 

The  last  approximation  was  applied  because  the  imaginary  part  of  the 


r:-^v?v:v  vivvi^ 


wavenumber  is  much  less  than  the  real  part  in  all  cases  of  interest. 
Considering  Eq.  (4.11),  one  can  see  that  even  if  the  object  inhomogeneities  are 
negligible  (k  =  k0),  the  difference  between  the  attenuation  coefficients  of  the 
object  and  medium  (c0  -  c(7))  may  make  the  imaginary  part  (and  thus  the 
magnitude)  of  the  object  function  large  enough  so  that  the  Born  approximation 
would  break  down.  Therefore,  we  expect  poor  reconstructions  for  large 
differences  in  the  attenuation  coefficient,  even  for  very  small  differences  in  the 
refractive  index  of  the  object  and  medium.  Perhaps  this  would  not  happen  in 
microwave  imaging  because  the  attenuation  constants  of  water-based  materials 
are  very  close.  It  should  be  noted  that  since  the  estimation  parameter  in  Eq. 
(4.11)  is  a  complex  function,  any  error  in  either  the  refractive  index  (k(7)/k0)  or 
the  attenuation  (c(7))  would  affect  the  estimated  value  of  the  other. 

In  all  the  simulation  results  in  later  sections,  we  will  refer  to  the  refractive 
index  as  the  ratio  of  k(7)  to  k0  (the  ratio  of  the  speed  of  the  wave  in  the 
medium  to  that  in  the  object),  rather  than  the  complex  function  nfr)  from  Eq. 
(4.9). 

4.2.1  A  Simplified  View  of  the  Problem 

In  chapter  10  we  showed  that  due  to  the  lack  of  an  inverse  Laplace 
transform  matching  the  domain  of  the  coverage,  the  current  reconstruction 
techniques  can  not  be  extended  to  the  case  of  an  attenuating  medium.  If  one  is 
not  dealing  with  high  attenuation  in  the  medium,  we  will  show  that  it  is 
possible  to  approximate  the  Laplace  transformation  by  a  Fourier  transfomation 
over  a  sub-domain  of  the  coverage.  We  will  consider  both  the  conventional 
technique  and  the  Synthetic  Aperture  technique. 

Consider  Eq.  (3.28)  derived  for  the  conventional  technique,  expressing  the 
relation  between  the  Fourier  transform  of  the  scattered  field  on  the  receiver 
line  and  the  object's  Laplace  Transform.  The  Laplace  transformation  can  be 
approximated  by  a  Fourier  transformation  only  if  the  imaginary  part  of 

A  -  kcoi/  =  (  a  ,  /j  -  kco  )  (4.12) 

is  negligible  (see  Eq.  (3.28)).  Since  a  is  real,  only  the  imaginary  part  of  J~kro 
has  to  be  small.  Let's  define  a  and  fn  as  fractions  of  k0 

o  =  a  kr,  ,  (4.13) 

C  =  b  kQ  (4. It) 

For  microwave  in  tissue,  the  upper  limit  on  b  is  about  0.07.  .?-kco  can  then  be 
approximated  by 
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3  ~  krn  -  J k  2  -  or"  -  kr 


1  -  a-  +  j  2  b  -  kc 


(4.15) 


Limiting  a  to  values  of  less  than  0.7,  approximately.  Eq.  (4.15)  can  still  be 
further  simplified 


-  a-  -  1 


3  ~  kco  -  ko  (  V  1  ~  a-  -  1  )  +  j  a'  £0  (4.16) 

This  restriction  on  a  is  not  required,  and  the  only  reason  we  have  applied  it  is 

to  make  the  mathematics  much  simpler.  The  imaginary  part  of  /?- kco  in  the 
equation  above  is  linearly  proportional  to  a2  or  in  fact  a2.  Therefore,  for  values 
of  a  smaller  than  an  upper  bound,  the  imaginary  part  of  3~ kco  is  negligible  and 
the  Laplace  transformation  in  Eq.  (3.28)  can  be  approximated  by  the  Fourier 
transformation.  This  point  has  been  mentioned  in  an  unpublished  work  of 
Mueller  [Mueller]. 

In  the  following  we  look  at  an  interesting  concept,  which  shows  that  the 
limitation  on  ci  as  suggested  in  the  last  paragraph  is  indeed  a  physical 
constraint  on  a  in  an  attenuating  medium.  As  we  mentioned  in  section  2.1,  the 

angular  spectrum  of  a  field  on  the  line  x  =  x,  propagates  to  line  x  =  x2 

according  to  the  following  relation 


A(k  x2)  =  A(k  X|)  eJ 


(4.17) 


Fig.  4.32  shows  a  plot  of  the  magnitude  of  the  exponential  factor  in  Eq.  (4.17) 
as  a  function  of  both  the  attenuation  of  the  medium  and  the  spatial  frequency 
of  the  plane  wave  in  the  y  direction.  We  have  plotted  the  magnitude  in  db  as 
a  function  of  the  dimensionless  parameter  7,  where  kx=7k0.  Thus,  for  7=0 
the  wave  is  traveling  directly  towards  the  receiver  line,  while  for  q~l  the  wave 
is  propagating  along  the  direction  of  the  receiver  line.  In  this  plot  the 
attenuation  factor  and  7  have  been  changed  from  0  to  5.0  db/cm  and  0  to  1.0, 
respectively.  It  is  clear  that  the  high  frequency  components  (larger  7)  are  more 
attenuated  than  those  components  at  lower  frequencies  (those  that  travel 
directly  towards  the  receiver  line).  This  means  that  as  the  field  propagates  in 
an  attenuating  medium,  it  loses  its  high  frequency  components. 

Remembering  the  Fourier  domain  coverage  of  diffraction  transmission 
tomography,  one  can  now  associate  this  phenomenon  with  a  degradation  in 
resolution.  This  point  is  illustrated  in  Fig.  4.33.  As  in  the  case  of  no 
attenuation,  the  inner  circle  corresponds  to  transmission  tomography,  while  the 
outer  ring  represents  the  data  measured  by  a  reflection  tomographic 
experiment.  The  difference  made  by  attenuation  is  shown  as  the  shaded  area. 
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Figure  t 32 

Energy  loss  of  plane  wave  components  resulting  from  an  angular  spectrum 
expansion  in  a  lossy  medium  as  a  function  of  directional  angle  and  at¬ 
tenuation  constant  of  the  medium.  7  is  a  measure  of  the  direction  of  pro¬ 
pagation  of  these  plane  waves.  7  =  0  corresponds  to  the  perpendicular 
direction  to  the  line  of  expansion,  and  7=1  corresponds  to  the  parallel 
direction  to  that  line. 


Figure  4.33 

Object’s  Fourier  domain  coverage  for  the  conventional  tomography  in 
transmission  and  reflection  modes,  in  a  lossy  medium.  The  lower  resolu¬ 
tion  in  the  conventional  transmission  tomography  is  due  to  the  smaller  ra¬ 
dius  of  the  coverage  circle  as  shown  here. 


In  this  region  the  attenuation  of  the  medium  reduces  the  amplitude  of  the 
plane  wave  components  below  a  minimum  tolerable  Signal  to  Noise  Ratio 
(SNR)  and  thus  makes  this  region  of  the  object's  Fourier  transform 
unmeasurable. 

In  the  Synthetic  Aperture  technique,  the  concept  of  propagation  of  the 
angular  spectrum  in  an  attenuating  medium  can  be  applied  again,  and  this 
results  in  upper  bounds  on  ky  and  J.  ky  is  the  spatial  frequency  associated 
with  the  transmitter,  and  3  is  the  spatial  frequency  associated  with  the 
receiver.  This  is  very  similar  to  the  effect  of  the  finite  size  of  the  transmitter 
and  receiver  arrays,  as  shown  in  Fig.  4.24.  This  point  is  illustrated  in  Fig. 
4.34.  The  coverage  in  the  reflection  mode  also  shrinks  as  shown  in  the  figure. 


4.2.2  Attenuation  in  Conventional  Diffraction  Tomography 

In  the  absence  of  a  reconstruction  procedure  specifically  designed  for  an 
object  in  a  lossy  medium  a  study  of  the  performance  of  the  present  algorithms 
in  the  presence  of  attenuation  is  in  order.  Since  the  quantitative  aspect  of  the 
reconstruction  is  of  importance  to  us,  we  have  compensated  for  the  loss  of 
energy  in  the  propagation  of  the  field  to  the  receiver  array.  In  the  following  we 
present  reconstructions  achieved  in  this  way.  based  on  the  projection  data 
calculated  using  the  programs  discussed  in  the  previous  sections. 

In  Figs.  4.35-47  we  have  shown  reconstructed  values  for  the  refractive 
index  and  the  attenuation  for  a  number  of  different  cylinders  and  media.  In 
each  case  the  first  plot  (a)  represents  the  difference  in  the  refractive  index  of 
the  cylinder  and  the  surrounding  medium  and  the  bottom  plot  (b)  shows  the 
reconstructed  difference  in  the  attenuation  constants.  In  the  first  set  of  Figs. 
4.35-40,  the  attenuation  constant  of  the  object  and  medium  are  taken  to  be 
equal.  Considering  these  figures,  we  conclude  that  higher  values  of  attenuation 
result  in  the  loss  of  higher  frequencies  in  the  reconstruction  and  result,  in  a  loss 
of  resolution.  The  non-zero  estimates  of  the  attenuation  constant  are  due  to 
the  error  in  the  phase  of  the  object  function. 

In  Figs.  4.41  through  4.44  we  have  considered  different  values  for  the 
attenuation  constant  of  the  object  am!  the  medium,  while  in  all  cases  we  have 
kept  the  refractive  index  of  the  object  and  the  medium  equal.  Again,  higher 
attenuation  of  the  medium  causes  the  loss  of  high  frequency  components  in  the 


reconstruction. 


Now,  let’s  examine  a  lossy  object  in  a  lossless  medium.  In  Figs.  4.45  and 
4.46.  we  have  shown  reconstructions  of  both  the  refractive  index  and  the 
attenuation  constant  of  cylinders  with  refractive  index  of  1.01  and  attenuation 
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Figure  *-3t 

When  medium  is  lossy,  larger  number  of  views  in  Synthetic  Aperture  tech¬ 
nique  can  increase  the  resolution  and  decrease  the  distortion  in  the  recon¬ 
structed  image  The  coverage  shown  corresponds  to  three  rotational  views 
of  object,  00°  apart. 
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Figure  4  35 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.01,  using  the  conventional  technique.  Both  the  object  and  the  medium 
are  lossless.  The  projections  are  sampled  at  64  points  over  a  sampling  in¬ 
terval  of  0.45  X.  The  receiver  line  is  at  a  distance  of  10  X  from  the 
cylinder,  (a)  Refractive  index,  (b)  Difference  in  attenuation  constant  of 
the  object  and  the  medium. 
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Figure  \  .55 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.01,  using  the  conventional  technique.  The  attenuation  constant  is  1.0 
db/cm  in  both  the  object  and  the  medium.  The  projections  are  sampled 
at  64  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  at  a 
distance  of  10  X  from  the  cylinder,  (a)  Refractive  index  (b)  Difference  in 
attenuation  constant  of  the  object  and  the  medium. 
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Figure  4.37 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.01,  using  the  conventional  technique.  The  attenuation  constant  is  3.0 
db/cm  in  both  the  object  and  the  medium.  The  projections  are  sampled 
at  64  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  at  a 
distance  of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference  in 
attenuation  constant  of  the  object  and  the  medium. 


<  a 


<■  b  ) 


Figure  4  38 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.03,  using  the  conventional  technique.  Both  the  object  and  the  medium 
are  lossless.  The  projections  are  sampled  at  64  points  over  a  sampling  in¬ 
terval  of  0.45  X.  The  receiver  line  is  at  a  distance  of  10  X  from  the 
cylinder,  (a)  Refractive  index,  (b)  Difference  in  attenuation  constant  of 
the  object  and  the  medium. 


I  f  ife 


15  5 


155  -ISO  A 


y  “CD 


stu 


-  mi  -^feg 


15  5^.lt  o 


Figure  4.40 

Reconstruction  of  a  cylinder  with  diameter  of  8  X  and  refractive  index  of 
1.03,  using  the  conventional  technique.  The  attenuation  constant  is  3.0 
db/cm  in  both  the  object  and  the  medium.  The  projections  are  sampled 
at  84  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  at  a 
distance  of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference  in 
attenuation  constant  of  the  object  and  the  medium. 
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Figure  -4.41 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.0,  using  the  conventional  technique.  The  attenuation  constant  of  the 
medium  and  the  object  are  0.1  db/cm  and  0.2  db/cm,  respectively.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X. 
The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder,  (a)  Refractive 
index,  (b)  Difference  in  attenuation  constant  of  the  object  and  the  medi- 
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Figure  4.42 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.0,  using  the  conventional  technique.  The  attenuation  constant  of  the 
medium  and  the  object  are  2.0  db/cm  and  2.1  db/cm,  respectively.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X. 
The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder,  (a)  Refractive 
index,  (b)  Difference  in  attenuation  constant  of  the  object  and  the  medi- 
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Figure  4.43 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.0,  using  the  conventional  technique.  The  attenuation  constant  of  the 
medium  and  the  object  are  0.1  db/cm  and  1.0  db/cm,  respectively.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X. 
The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder,  (a)  Refractive 
index,  (b)  Difference  in  attenuation  constant  of  the  object  and  the  medi¬ 
um. 


Figure  I  II 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.0,  using  the  conventional  technique.  The  attenuation  constant  of  the 
medium  and  the  object  are  2.0  db/cm  and  3.0  db/cm,  respectively.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X. 
The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder,  (a)  Refractive 
index,  (b)  Difference  in  attenuation  constant  of  the  object  and  the  medi¬ 
um. 


constants  of  0.5  and  3.0  db/cm,  respectively.  It  is  only  for  the  ca^c  of  3.0 
db/cm  attenuation  that  the  reconstruction  is  really  distorted  quantitatively, 
while  for  the  value  of  0.5  db/cm  attenuation  we  get  a  reasonably  accurate 
reconstruction.  Since  in  most  practical  cases  of  interest  the  difference  between 
the  attenuation  constants  of  the  object  and  the  medium  is  not  larger  than  10 
db/cm,  perhaps  the  difference  in  the  attenuation  constants  is  not  a  limiting 
factor  for  the  reconstruction  algorithms  in  a  lossy  medium. 

In  the  final  reconstruction  (Fig.  4.17)  we  have  tried  to  consider  a  case  close 
to  reality,  taking  all  the  parameters  of  the  object  and  medium  to  be  different. 
The  refractive  index  of  the  object  has  been  set  to  1.02  and  the  difference  in 
attenuation  to  0.5  db/cm.  The  estimation  error  in  the  phase  of  the  object 
function  causes  the  estimate  of  the  refractive  index  to  become  smaller  than  the 
actual  value,  while  the  difference  in  attenuation  is  overestimated. 

The  simulation  results  presented  in  this  section  show  that  the  only  effect 
of  the  medium  attenuation  in  the  conventional  approach  is  some  loss  of 
resolution. 

4.2.3  Attenuation  in  Synthetic  Aperture  Diffraction  Tomography 

As  discussed  in  chapter  10,  one  can  not  extend  the  Synthetic  Aperture 
technique  to  the  case  of  an  attenuating  medium.  Therefore,  in  the  absence  of  a 
reconstruction  procedure  specifically  designed  for  a  lossy  medium  we  will 
examine  the  performance  of  the  present  algorithm  in  the  presence  of 
attenuation.  As  in  the  conventional  tomography,  the  quantitative  aspect  of 
the  reconstruction  is  of  importance  to  us,  thus  we  should  compensate  for  the 
loss  of  energy  in  the  propagation  of  the  field  to  the  receiver  array.  However, 
the  compensation  is  quite  different  from  that  of  the  conventional  technique.  In 
the  conventional  technique,  the  attenuation  is  compensated  for  by  multiplying 
the  each  projection  data  by  a  constant  value,  while  in  tne  bynthetic  Aperture 
technique  the  multiplication  factor  is  dependent  on  the  transducers  positions. 
The  scattered  field  at  each  position  of  the  transmitter  and  receiver  is  amplified 
by  the  amount  of  attenuation  the  incident  field  has  suffered  traveling  between 
the  two  transducers.  This  means  that  the  field  at  the  end  of  the  arrays  is 
weighted  higher  relative  to  the  one  at  the  renter.  Of  course,  in  the  presence  of 
noise,  the  noise  will  be  amplified  as  well.  In  the  following,  we  present 
reconstructions  achieved  in  this  way,  based  on  the  projection  data  calculated 
using  the  programs  discussed  in  chapter  2. 
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Figure  4.45 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.01,  using  the  conventional  technique.  The  medium  is  lossless  while  the 
attenuation  constant  of  the  object  is  0.5  db/cm.  The  projections  are  sam¬ 
pled  at  64  points  over  a  sampling  interval  of  0.45  X.  The  receiver  line  is  at 
a  distance  of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference 
in  attenuation  constant  of  the  object  and  the  medium. 
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Figure  4  16 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.01,  using  the  conventional  technique.  The  medium  is  lossless  while  the 
attenuation  constant  of  the  object  is  3.0  db/cm.  The  projections  are  sam¬ 
pled  at  64  points  over  a  sampling  interval  of  0  45  X.  The  receiver  line  is  at 
a  distance  of  10  X  from  the  cylinder,  (a)  Refractive  index,  (b)  Difference 
in  attenuation  constant  of  the  object  and  the  medium. 
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Figure  4  17 

Reconstruction  of  a  cylinder  with  diameter  of  6  X  and  refractive  index  of 
1.02,  using  the  conventional  technique.  The  attenuation  constants  of  the 
medium  and  the  object  are  3.0  db/cm  and  3.5  db/cm,  respectively.  The 
projections  are  sampled  at  64  points  over  a  sampling  interval  of  0.45  X. 
The  receiver  line  is  at  a  distance  of  10  X  from  the  cylinder,  (a)  Refractive 
index,  (b)  Difference  in  attenuation  constant  of  the  object  and  the  medi¬ 
um. 
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In  Fig.  4.48a,b  we  have  shown  the  magnitude  of  the  scattered  field  from  a 
cylinder  with  a  refractive  index  of  1.03  and  a  diameter  of  6  wavelengths  and  its 
Fourier  transform.  The  attenuation  constants  of  the  medium  and  the  object 
were  set  to  3.0  db/cm.  As  one  expects  intuitively,  the  field  goes  to  zero  at  the 
end  of  the  receiver  array  and  if  no  attenuation  compensation  is  performed,  the 
length  of  the  arrays  would  not  affect  the  reconstruction.  From  Fig.  4.48b,  the 
attenuation  of  high  frequency  components  is  obvious.  The  reconstructed 
object’s  Fourier  transform  is  shown  in  Fig.  4.48c.  It  is  quite  clear  that  those 
values  in  the  diagonal  directions  are  attenuated  and  will  cause  some  distortion 
in  those  directions.  However,  if  the  attenuation  is  compensated  for,  this  effect 
will  be  reduced.  In  Fig.  4.49  we  have  shown  the  magnitude  of  the 
reconstructed  a  refractive  index  of  a  cylinder  with  refractive  index  of  1.03  and 
diameter  of  6  wavelengths  when  the  attenuation  is  compensated  for.  Tim 
attenuation  constant  of  the  medium  and  the  object  are  both  taken  to  be  3.0 
db/cm.  In  Fig.  4.50  we  have  shown  the  reconstruction  when  no  attenuation 
compensation  has  been  performed.  Of  course  the  reconstruction  is  qualitatively 
wrong  and  some  distortion  is  present  along  the  diagonal  directions.  When  four 
views,  45  degrees  apart,  are  used  the  distortion  is  removed  as  shown  in  Fig. 
4.51.  In  Figs.  4.50  and  4.51  the  smoothness  of  the  reconstruction  is  due  to  the 
loss  of  high  frequency  components  in  the  projection  due  to  the  medium 
attenuation. 

In  Fig.  4.52,  the  reconstructed  difference  between  attenuation  constant  of 
the  object  and  medium  is  presented.  The  attenuation  constants  of  the  medium 
and  the  object  were  set  to  2.0  and  3.0  db/cm,  respectively.  The  cylinder  was  6 
wavelengths  in  diameter  and  the  refractive  index  of  the  object  and  medium 
were  taken  to  be  equal.  Other  than  a  small  amount  of  distortion  in  the 
diagonal  directions,  the  quality  of  the  reconstruction  is  satisfactory. 

Finally,  in  Fig.  4.53  we  have  shown  the  reconstruction  results  for  the 
object  considered  in  Fig.  4.47.  The  result  obtained  under  the  Synthetic 
Aperture  technique  (Fig.  4.53)  is  quite  similar  to  the  one  obtained  using  the 
conventional  technique  (Fig.  4.47).  In  both  of  these  reconstructions  the 
attenuation  has  been  compensated  for. 
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Figure  4.48 

The  object  is  a  lossy  cylinder  with  diameter  and  refractive  index  of  8  X 
and  1.03,  respectively.  The  attenuation  constants  of  the  medium  and  the 
object  are  both  3.0  db/cm.  There  are  128  elements  in  the  transmitter  and 
the  receiver  arrays,  with  the  inter-element  distance  of  0.45  X.  (a)  Magni¬ 
tude  of  the  scattered  field,  as  a  function  of  the  transmitter  and  receiver 
positions,  (b)  Magnitude  of  the  Fourier  transform  of  the  scattered  field, 
with  respect  to  the  transmitter  and  receiver  positions,  (c)  Magnitude  of 
the  reconstructed  object’s  Fourier  transform. 
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Figure  150 

Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with  diame¬ 
ter  of  8  X  and  refractive  index  of  1.03,  using  the  Synthetic  Aperture  tech¬ 
nique.  The  attenuation  constant  of  the  object  and  the  medium  are  both 
3.0  db/cm.  The  transmitter  and  the  receiver  arrays,  positioned  20  X 
apart,  nave  128  elements  with  the  inter-element  distance  of  0.45  X.  The 
attenuation  of  the  medium  has  not  been  compensated  for. 
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Figure  4.oI 

Magnitude  of  the  reconstructed  refractive  index  of  a  cylinder  with  diame¬ 
ter  of  6  X  and  refractive  index  of  i.03,  using  the  Synthetic  Aperture  tech¬ 
nique.  The  attenuation  constant  of  the  object  and  the  medium  are  both 
3.0  db/cm.  The  transmitter  and  the  receiver  arrays,  positioned  20  X 
apart,  nave  128  elements  with  the  inter-element  distance  of  0.45  X.  The 
attenuation  of  the  medium  has  not  been  compensated  for,  but  four  rota¬ 
tional  views  of  the  object  has  been  used. 
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Figure  4.52 

Magnitude  of  the  reconstructed  attenuation  constant  of  a  cylinder  with  di¬ 
ameter  of  6  X  and  refractive  index  of  1.0,  using  the  Synthetic  Aperture 
technique.  The  attenuation  constants  of  the  object  and  medium  are  3.0 
db/cm  and  2.0  db/cm,  respectively.  The  transmitter  and  the  receiver  ar¬ 
rays,  positioned  20  X  apart,  have  128  elements  with  the  inter-element  dis¬ 
tance  of  0.45  X.  The  attenuation  of  the  medium  has  been  compensated 
for. 


List  of  References 


[Mueller]  R.  K.  Mueller  and  M.  Kaveh,  “Diffraction  tomography,”  Unpub¬ 
lished  work. 

[Nahamoo82]  D.  Nahamoo  and  A.  C.  Kak,  "Ultrasonic  diffraction  imaging,” 
Purdue  University,  School  of  Electrical  Engineering,  Technical 
Report,  TR-EE-82-20,  1982. 

[Nahamoo84]  D.  Nahamoo,  S.  X.  Pan,  and  A.  C.  Kak,  “Synthetic  aperture 
diffraction  tomography  and  its  interpolation-free  computer  im¬ 
plementation,”  IEEE  Transactions  on  Sonics  and  Ultrasonics, 
Vol.  SU-31,  No.  4,  pp.  218-229,  July  1984. 


I 


CHAPTER  5 
SUMMARY 


To  give  the  reader  a  brief  overview  of  the  work  presented  in  the  second 
part  of  this  thesis,  we  have  devoted  this  chapter  to  the  summary  of  the  topics 
and  results  already  discussed  in  chapters  2,  3,  and  4. 

The  forward  problem  in  diffraction  imaging,  by  which  we  mean  calculation 
cf  mattered  fields  from  the  inhomogeneities  present  in  an  object,  was  addressed 
in  chapter  2.  Two  distinct  types  of  objects  were  considered:  single-component 
objects  and  multi-component  objects.  The  object  components  were  assumed  to 
be  cylinders,  because  in  the  formulation  of  the  problem  for  non-cylindrical 
objects  one  is  faced  with  severe  mathematical  difficulties. 

For  incident  plane  waves,  we  used  the  well-known  mathematical 
expressions  for  the  scattered  fields  derived  in  most  of  the  classical  books  on 
scattering  theory.  For  the  case  of  an  incident  field  generated  by  a  point  source, 
we  derived  similar  expressions  for  the  exact  scattered  fields  using  the  series 
expansion  of  a  Hankel  function  in  cylindrical  coordinates.  A  computer  program 
was  developed  to  simulate  the  process  by  utilizing  the  mathematical 
formulations.  Some  parts  of  the  code  was  vectorized  to  speed  up  the  floating 
point  calculations  in  the  array  processor  FPS-AP120B.  By  using  the  developed 
computer  program,  we  are  now  capable  of  generating  the  exact  scattered  fields 
from  a  cylinder  with  arbitrary  size,  refractive  index,  and  attenuation  constant 
embedded  in  a  medium  with  an  arbitrary  attenuation  constant,  illuminated  by 
either  a  plane  wave  or  a  point  source. 

For  calculating  the  scattered  fields  from  multi-component  objects,  a  major 
source  of  difficulty  is  dealing  with  the  interaction  between  the  various 
components.  Depending  on  the  interaction  between  the  components,  the  total 
scattered  field  may  or  may  not  bear  any  resemblance  to  the  simple  sum  of  the 
scattered  fields  for  each  of  the  components,  assuming  the  others  to  be  absent. 
In  chapter  2,  we  presented  a  new  computational  procedure  for  calculating  the 
inter-component  interaction.  We  have  based  our  technique  on  a  basic  theory  in 
which  the  scattered  field  is  expressed  as  an  infinite  summation  of  different 
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order  terms.  For  the  first  order  scattering,  each  component  interacts 
independently  with  the  incident  illumination,  being  oblivious  of  the  existence  of 
the  others.  Second  order  scattering  terms  are  generated  when  each  component 
interacts  with  the  first  order  fields  sent  in  its  direction  by  the  other 
components,  and  so  on.  Second  order  fields  are  calculated  by  expanding  the 
first  order  fields  from  each  component  as  superposition  of  plane  waves  and  then 
analyzing  the  interaction  of  each  of  these  plane  waves  with  the  other  object 
components.  Higher  order  fields  can  be  generated  similarly. 

This  technique  was  coded  in  a  vectorized  fashion  to  benefit  from  the  speed 
of  available  array  processor.  Although  our  program  only  calculates  first  and 
second  order  scattering  terms,  extension  of  the  code  for  calculation  of  higher 
order  terms  is  trivial.  Using  this  procedure  we  showed  by  computer 
simulations  that  even  when  each  component  of  a  two-component  object  is 
weakly  scattering,  the  multiple  scattering  effects  become  important  when  the 
components  are  blocking  each  other.  We  further  showed  that  when  strongly 
scattering  components  that  are  large  compared  to  a  component  are  not 
blocking  each  other,  the  scattering  effects  can  be  ignored.  Both  these 
conclusions  agree  with  intuitive  reasoning. 

In  chapter  3.  we  investigated  the  inverse  scattering  problem.  In  general 
inverse  scattering  is  a  much  more  complicated  problem  than  the  forward 
problem.  Inverse  scattering  techniques  are  usually  accompanied  by  severe 
restrictions  on  the  different  parameters  of  the  system  under  consideration.  We 
have  considered  those  inverse  scattering  techniques  which  are  based  on  the 
Fourier  Diffraction  Projection  Theorem.  This  theorem  related  the  Fourier 
transform  of  an  object  with  the  Fourier  transform  of  the  scattered  field 
measured  on  a  straight  receiver  line  when  the  object  is  illuminated  by  a  plane 
wave.  It  should  be  noted  that  the  Fourier  Diffraction  Projection  Theorem  is 
valid  only  when  the  object  inhomogeneities  are  weakly  scattering,  or  in  fact 
when  either  the  two  well  known  Born  or  Rytov  approximations  hold. 

The  two  inverse  reconstruction  techniques  considered  are:  the  conventional 
technique  and  the  Synthetic  Aperture  technique.  In  the  conventional  technique 
the  object  is  illuminated  by  a  plane  wave  and  the  scattered  field  is  measured  on 
a  receiver  line  which  is  perpendicular  to  the  direction  of  the  propagation  of  the 
incident  field.  Many  rotational  views  of  the  object  provides  the  necessary  data 
to  fill  up  the  Fourier  space  of  the  object  inside  a  circle  of  radius  v^kr,  (k0  is  the 
wavenumber).  The  object's  cross  section  is  then  recovered  by  an  inverse 
Fourier  transformation 


The  Synthetic  Aperture  technique  is  a  more  general  technique  which 
permits  the  use  of  any  type  of  insonification,  as  opposed  to  only  plane  waves 
and  also  requires  only  two  rotational  view's  of  the  object.  However,  for  each 
rotational  position,  the  projection  data  is  collected  by  independent  traversal 
movements  of  the  transmitting  and  the  receiving  transducers,  over  two  parallel 
lines  on  both  sides  of  the  object.  In  this  technique,  the  coverage  in  the  Fourier 
domain  of  the  object  is  inside  a  circle  of  radius  v/2k0,  the  same  as  that  of  the 
conventional  approach. 

We  have  modified  the  theory  behind  the  conventional  reconstruction 
technique  to  include  the  case  of  lossy  objects  and  media  as  well.  We  have 
shown  that  in  a  lossy  medium  one  would  recover  the  two  dimensional  Laplace 
transform  of  the  object  over  a  surface  and  moreover  the  inverse  Laplace 
transform  techniques  can  not  be  used  to  recover  the  object’s  cross  section. 
Using  the  Synthetic  Aperture  technique  in  a  lossy  medium,  the  object's 
Fourier-Laplace  transform  is  recovered  over  a  surface  in  the  Fourier-Laplace 
domain  of  the  object.  By  Fourier-Laplace  transformation  we  mean  Fourier 
transformation  in  one  direction  and  Laplace  transformation  in  the  other 
direction.  Still  no  inverse  transform  technique  can  be  used  for  the  inversion 
process.  A  possible  way  to  avoid  the  inverse  Fourier-Laplace  transformation  is 
to  impose  a  constraint  on  the  processed  data  to  convert  the  transform  to  a 
simple  two-dimensional  Fourier  transform.  However,  the  experimental 
procedure  has  to  be  repeated  for  many  different  rotational  views  of  the  object. 
Of  course  this  requires  enormous  amount  of  data  to  be  collected  and  may  be 
quite  restrictive  in  some  applications. 

In  chapter  3,  it  has  also  been  proven  that  if  the  Synthetic  Aperture 
technique  is  used  in  the  reflection  mode  as  well  as  the  transmission  mode,  one 
rotational  position  of  the  object  is  enough  for  recovering  the  object’s  Fourier 
transform.  In  addition,  the  coverage  would  be  inside  a  circle  of  radius  2k0 
rather  than  v/2k0  for  the  conventional  and  the  Synthetic  Aperture  techniques. 
Further  study  of  this  technique  still  remains  to  be  examined  by  simulation.il 
means. 

Chapter  1  is  devoted  to  a  simulation  study  of  multiple  scattering  and 
attenuation  phenomena  in  diffraction  imaging.  The  programs  developed  in 
chapter  2  Ur  the  calculation  of  scattered  fields  together  with  the  reconstruction 
schemes  presented  in  chapter  3  were  used  to  generate  the  computer  simulation 
results  in  this  chapter. 


All  diffraction  tomography  algorithms  available  today,  including  the 
conventional  and  the  Synthetic  Aperture  techniques,  assume  the  first  order 
Bom  or  Rytov  approximations;  which  implies  that  they  are  only  good  for 
reconstructing  weakly  scattering  objects.  In  a  multi-component  object  this 
translates  into  ignoring  the  interaction  between  the  components.  However  in 
most  cases  of  practical  interest  the  object  would  be  composed  of  several 
distinct  components.  One  of  our  objectives  in  chapter  4  has  been  to  study  the 
effect  of  the  interaction  between  object  components  on  the  reconstructed 
object,  in  both  the  conventional  and  the  Synthetic  Aperture  techniques. 

We  generated  two  series  of  reconstructions  of  a  multi-cylinder  object  using 
the  conventional  reconstruction  technique.  In  the  first  series,  the  projections 
were  generated  by  ignoring  the  second  and  higher  order  scattered  fields.  In  the 
other  the  second  order  scattered  fields  were  included  in  the  projections.  A 
comparison  of  these  two  series  of  simulations  showed  that  the  scattering 
interaction  between  the  components  cause  some  distortion  along  the  line 
joining  the  components,  the  distortion  increasing  with  larger  values  of  the 
refractive  indices  of  the  components.  The  distortion  also  becomes  more 
noticeable  as  the  number  of  components  increases. 

We  examined  three  different  approaches  for  reducing  the  distortions  in  the 
conventional  technique.  First  we  made  the  receiver  array  to  be  oblique  to  the 
direction  of  the  incident  field.  In  this  way  the  blocking  effect,  which  is  quite 
severe  when  the  propagation  vector  and  the  line  joining  the  object  components 
are  collincar,  is  reduced.  Simulation  results  showed  that  although  this 
arrangement  reduces  the  distortion  in  the  region  between  the  object 
components,  the  price  paid  is  an  increase  in  the  general  distortion  all  through 
the  object.  This  approach  seems  appropriate  for  two-component  objects. 

In  the  second  approach  intended  for  more  complex  objects,  the  object  is 
illuminated  with  several  non-collinear  plane  waves  for  each  projection,  for  each 
rotational  position  of  the  receive  array.  Fach  plane  wave  generates  a  complete 
frequency  domain  coverage,  these  are  averaged  prior  to  Fourier  inversion  for 
object  reconstruction.  The  rationale  for  this  method  is  that  the  error  caused 
by  blocking  vary  for  the  three  different  illuminating  waves.  The  simulation 
results  were  similar  to  that  of  the  first  approach. 

We  also  considered  another  approach  which  is  similar  in  spirit  to  the 
Synthetic  Aperture  technique  For  each  position  of  the  receiver  line  the  object 
is  illuminated  by  incident  plane  waves  from  several  directions.  However,  unlike 
the  second  approach,  no  averaging  is  performed  in  the  frequency  domain.  The 
object  is  now  rotated  by  a  large  angle  to  yield  a  similar  coverage  over  areas  left 
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blank  before;  this  procedure  is  repeated  until  no  holes  are  left  in  the  frequ 
domain.  The  rationale  here  is  that  with  each  plane  wave  at  a  difff 
orientation  with  respect  to  the  receiver  line,  the  severity  of  the  blocking  e 
caused  by  multiple  scattering  should  be  reduced.  Simulation  results  show 
this  approach  does  result  in  better  reconstructions  of  multi-component  obje< 

The  superiority  of  the  third  approach  to  the  others  persuaded  us  that 
performance  of  the  Synthetic  Aperture  technique  in  the  presence  of  m 
component  interactions  should  be  superior  to  that  of  the  conventh 
technique.  Simulation  results  proved  that  our  expectation  was  in  fact  t 
There  was  almost  no  distortion  between  the  components  as  was  the  case  in 
reconstructions  achieved  with  the  conventional  technique.  The  reason  is  t 
the  blocking  effect  is  minimized  in  this  technique.  Another  interesting  p( 
was  that  the  components  which  were  blocked  became  more  distorted  compa 
to  the  ones  at  closer  distances  to  the  transmitter  array.  Our  general  conclus 
was  that  the  performance  of  the  Synthetic  Aperture  technique  in  the  prese 
of  multiple  scattering  is  far  superior  than  that  of  the  conventional  technique. 

In  the  second  half  of  the  chapter  4,  we  discussed  the  effect  of 
attenuation  on  the  object  reconstructions  in  diffraction  tomography.  Althoi 
in  ultrasonic  imaging  the  role  of  attenuation  is  minor,  in  microwave  imagin' 
can  not  be  ignored,  because  microwaves  at  1  GHz  undergo  3  db  of  attenuat 
per  centimeter  in  water.  In  the  lack  of  any  reconstruction  technique  specifi 
designed  for  lossy  media,  we  studied  the  effect  of  medium  attenuation  on 
calculated  projections  and  thereby  on  the  reconstructions  themselves. 

We  proved  that  the  medium  attenuation  suppresses  the  high  frequer 
components  present  in  the  projections.  In  the  conventional  technique  this  le: 
to  a  loss  of  resolution  in  the  final  object's  reconstruction.  This  point  is  back 
by  our  computer  simulation  results.  However  in  the  Synthetic  Aperti 
technique,  medium  attenuation  causes  a  loss  of  coverage  in  the  object’s  Four 
domain  in  a  manner  that  the  lost  region  can  be  recovered  only  by  consider; 
more  rotational  views  of  the  object.  Heuristic  compensation  of  medii 
attenuation  in  both  techniques  results  is  acceptable  quantitative  obj 
reconstructions.  However  such  compensation  techniques  may  not  be  advisal 
in  the  presence  of  noise  in  the  collected  projections. 

Finally  simultaneous  reconstruction  of  the  refractive  index  and  t 
attenuation  constant  of  an  object  was  discussed  in  chapter  4.  It  was  sho1 
that  large  refractive  index  of  an  object  could  result  in  an  erroner 
reconstruction  of  the  attenuation  constant,  while  typical  small  differem 
between  the  attenuation  constants  of  the  medium  and  the  object  would  r 
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affect  the  reconstructed  refractive  index  of  the  object.  As  long  as  the  noise  level 
is  under  certain  threshold  so  that  the  heuristic  compensations  for  the  medium 
attenuation  can  be  applied,  the  attenuation  of  the  medium  and  object  will  not 
be  a  limiting  factor. 
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APPENDIX 


The  aim  in  this  appendix  is  to  solve  the  following  set  of  equations  for  u;  in 
terms  of  ur,  v,  k0,  and  (Q 

(li+k,)2  +  (v+k,)2  =k  J  ,  (1) 

k„2  +  k,2  =  kto2  ,  (2) 

<1  -  u,  +  j  Uj  (3) 

k„  =  k„  +  j  <„  ■  (■») 

Substituting  Eq.  (2)  in  (1),  we  have: 

u2  +  v2  +  2vky  =  ~2ukx  .  (5) 

Squaring  both  sides  of  Eq.  (5)  and  making  use  of  Eq.  (2),  we  get 

u4  +  (v2  +  2vky)2  +  2u2(v2  +  2vky)  -  4u2(kco2-ksyby2)  =0  (6) 

Real  and  imaginary  part  of  Eq.  (6)  define  the  following  two  equations: 

(ur2-u2)2  -  4ur2u2  +  (v2  +  2vky)2  +  2(ur2~u2)  (v2  +  2vky) 

-  4(ur2-u;2)  (k02-t02~ky2)  +  16urUjk0t0  =  0  (7) 

uru,K2"ui2)  +  urUi(v2  +  2vky)  -  2(ur2-Uj2)k0e0 

~2urui(jc02-c02-ky2)  =  0  (8) 

Substitution  of  the  expression  for  ky  from  Eq.  (8)  into  Eq.  (7),  provides  us  with 

the  following  fifth  order  equation  for  U;  in  terms  of  ur,  v,  kc,  and  fc. 

76ui5  +  lb"*  +  I4"?  +  73ui2  +  TjV  +  7i  =  0  (9) 

where 

76  "  .  (Wa) 


72  =  -ur(ur4  +  v4  -  4(k02-«o2v2  +  2urV  ) 


(lOe 


*71  =  4ur2k0f0(ur2  +  v2)  (10f 

t 

\  second  order  equation  for  ky  is  achieved  by  rearranging  the  terms  in  Eq 


ky2  +  vky  +  7(Ur2-Ui2  +  V2)  -  (k02~f02)  -  -■-“'-koto  =  0  • 

i  UfU; 
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